diff --git a/NumLib/Extrapolation/ExtrapolatableElement.h b/NumLib/Extrapolation/ExtrapolatableElement.h
new file mode 100644
index 0000000000000000000000000000000000000000..6a666abe8cb4aa3555a09aace876ec725d5a2710
--- /dev/null
+++ b/NumLib/Extrapolation/ExtrapolatableElement.h
@@ -0,0 +1,35 @@
+/**
+ * \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
+
+#include <Eigen/Core>
+
+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::RowVectorXd> 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..432fab4bc1097728611580cc5337001ebe60a327
--- /dev/null
+++ b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
@@ -0,0 +1,134 @@
+/**
+ * \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).
+ *
+ * This is an interface class; suitable implementations have to be provided for
+ * concrete collections of extrapolatable elements.
+ */
+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::RowVectorXd> getShapeMatrix(
+ std::size_t const id, unsigned const integration_point) const = 0;
+
+ /*! Returns integration point values of some property of a specific element.
+ *
+ * \param id ID of the element of which the property is queried.
+ * \param cache Can be used to compute a property on-the-fly.
+ *
+ * \remark
+ * \parblock
+ * Since this method returns a reference, integration point values can not
+ * be computed locally and just returned.
+ *
+ * Instead, 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
+ */
+ 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 ExtrapolatableElement.");
+
+ /*! A method providing integration point values of some property.
+ *
+ * The method must return reference to a vector containing the integration
+ * point values.
+ *
+ * For further information about the \c cache parameter see
+ * ExtrapolatableElementCollection::getIntegrationPointValues().
+ */
+ 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::RowVectorXd> 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..cee424b32b95a49bb7d4f98330023660decbd3c3 100644
--- a/NumLib/Extrapolation/Extrapolator.h
+++ b/NumLib/Extrapolation/Extrapolator.h
@@ -14,64 +14,19 @@
#include <Eigen/Eigen>
#include "NumLib/NumericsConfig.h"
+#include "ExtrapolatableElementCollection.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 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 +35,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-impl.h
deleted file mode 100644
index 5e03eece5e3ba05e9c49ec5d82bf2ad31b270eba..0000000000000000000000000000000000000000
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
+++ /dev/null
@@ -1,153 +0,0 @@
-/**
- * \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
- *
- */
-
-#include <functional>
-
-#include <logog/include/logog.hpp>
-#include <Eigen/Core>
-
-#include "MathLib/LinAlg/LinAlg.h"
-#include "NumLib/Assembler/SerialExecutor.h"
-#include "MathLib/LinAlg/MatrixVectorTraits.h"
-#include "NumLib/Function/Interpolation.h"
-#include "LocalLinearLeastSquaresExtrapolator.h"
-
-namespace NumLib
-{
-
-template<typename PropertyTag, typename LocalAssembler>
-void
-LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
-extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
-{
- _nodal_values.setZero();
-
- // counts the writes to each nodal value, i.e., the summands in order to
- // compute the average afterwards
- auto counts =
- MathLib::MatrixVectorTraits<GlobalVector>::newInstance(_nodal_values);
- counts->setZero(); // TODO BLAS?
-
- using Self = LocalLinearLeastSquaresExtrapolator<
- PropertyTag, LocalAssembler>;
-
- NumLib::SerialExecutor::executeMemberDereferenced(
- *this, &Self::extrapolateElement, local_assemblers, property, *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)
-{
- assert(static_cast<std::size_t>(_residuals.size()) == local_assemblers.size());
-
- using Self = LocalLinearLeastSquaresExtrapolator<
- PropertyTag, LocalAssembler>;
-
- NumLib::SerialExecutor::executeMemberDereferenced(
- *this, &Self::calculateResiudalElement, local_assemblers, property);
-}
-
-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)
-{
- auto const& integration_point_values = loc_asm.getIntegrationPointValues(
- property, _integration_point_values_cache);
-
- // number of nodes in the element
- const auto nn = loc_asm.getShapeMatrix(0).cols();
- // number of integration points in the element
- const auto ni = integration_point_values.size();
-
- assert(ni >= static_cast<decltype(ni)>(nn) &&
- "Least squares is not possible if there are more nodes than"
- "integration points.");
-
- 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);
- assert(shp_mat.cols() == nn);
-
- // copy shape matrix to extrapolation matrix columnwise
- N.block(int_pt, 0, 1, nn) = shp_mat;
- }
-
- // 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());
-
- // TODO
- // optimization: Store decomposition of N*N^T or N^T for reuse?
- // cf. http://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html
- // cf. http://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html
- // cf. http://eigen.tuxfamily.org/dox/classEigen_1_1LDLT.html
- // cf. https://eigen.tuxfamily.org/dox/classEigen_1_1HouseholderQR.html
- // Extrapolate several values at once?
-
- // do the least squares computation using QR decomposition.
- Eigen::VectorXd tmp = N.householderQr().solve(integration_point_values_vec);
-
- // option: do the least squares computation using LDLT decomposition.
- // Eigen::VectorXd tmp =
- // (N*N.transpose()).ldlt().solve(N*integration_point_values_vec);
-
- // 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?
- 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)
-{
- auto const& gp_vals = loc_asm.getIntegrationPointValues(
- property, _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;
-
- // 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) {
- // TODO PETSc negative indices?
- nodal_vals_element[i] = _nodal_values[global_indices[i]];
- }
-
- double residual = 0.0;
- double gp_val_extrapol = 0.0;
-
- for (unsigned gp=0; gp<ni; ++gp)
- {
- NumLib::shapeFunctionInterpolate(
- nodal_vals_element, loc_asm.getShapeMatrix(gp), gp_val_extrapol);
- auto const& ax_m_b = gp_val_extrapol - gp_vals[gp];
- residual += ax_m_b * ax_m_b;
- }
-
- _residuals.set(element_index, std::sqrt(residual / ni));
-}
-
-} // namespace ProcessLib
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8979905afd346c007e4892d2dc8d35073a39dd42
--- /dev/null
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
@@ -0,0 +1,173 @@
+/**
+ * \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
+ *
+ */
+
+#include "LocalLinearLeastSquaresExtrapolator.h"
+
+#include <functional>
+
+#include <Eigen/Core>
+#include <logog/include/logog.hpp>
+
+#include "MathLib/LinAlg/LinAlg.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "NumLib/Assembler/SerialExecutor.h"
+#include "NumLib/Function/Interpolation.h"
+#include "ExtrapolatableElementCollection.h"
+
+namespace NumLib
+{
+LocalLinearLeastSquaresExtrapolator::LocalLinearLeastSquaresExtrapolator(
+ NumLib::LocalToGlobalIndexMap const& dof_table)
+ : _nodal_values(NumLib::GlobalVectorProvider::provider.getVector(
+ MathLib::MatrixSpecifications(dof_table.dofSizeWithoutGhosts(),
+ dof_table.dofSizeWithoutGhosts(),
+ &dof_table.getGhostIndices(),
+ nullptr)))
+#ifndef USE_PETSC
+ , _residuals(dof_table.size())
+#else
+ , _residuals(dof_table.size(), false)
+#endif
+ , _local_to_global(dof_table)
+{
+ /* Note in case the following assertion fails:
+ * If you copied the extrapolation code, for your processes from
+ * somewhere, note that the code from the groundwater flow process might
+ * not suit your needs: It is a special case and is therefore most
+ * likely too simplistic. You better adapt the extrapolation code from
+ * some more advanced process, like the TES process.
+ */
+ assert(dof_table.getNumberOfComponents() == 1 &&
+ "The d.o.f. table passed must be for one variable that has "
+ "only one component!");
+}
+
+void LocalLinearLeastSquaresExtrapolator::extrapolate(
+ ExtrapolatableElementCollection const& extrapolatables)
+{
+ _nodal_values.setZero();
+
+ // counts the writes to each nodal value, i.e., the summands in order to
+ // compute the average afterwards
+ auto counts =
+ MathLib::MatrixVectorTraits<GlobalVector>::newInstance(_nodal_values);
+ counts->setZero(); // TODO BLAS?
+
+ 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);
+}
+
+void LocalLinearLeastSquaresExtrapolator::calculateResiduals(
+ ExtrapolatableElementCollection const& extrapolatables)
+{
+ assert(static_cast<std::size_t>(_residuals.size()) ==
+ extrapolatables.size());
+
+ auto const size = extrapolatables.size();
+ for (std::size_t i=0; i<size; ++i) {
+ calculateResidualElement(i, extrapolatables);
+ }
+}
+
+void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
+ std::size_t const element_index,
+ ExtrapolatableElementCollection const& extrapolatables,
+ GlobalVector& counts)
+{
+ auto const& integration_point_values =
+ extrapolatables.getIntegrationPointValues(
+ element_index, _integration_point_values_cache);
+
+ // number of nodes in the element
+ const auto nn = extrapolatables.getShapeMatrix(element_index, 0).cols();
+ // number of integration points in the element
+ const auto ni = integration_point_values.size();
+
+ assert(ni >= static_cast<decltype(ni)>(nn) &&
+ "Least squares is not possible if there are more nodes than"
+ "integration points.");
+
+ 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 =
+ extrapolatables.getShapeMatrix(element_index, int_pt);
+ assert(shp_mat.cols() == nn);
+
+ // copy shape matrix to extrapolation matrix columnwise
+ N.block(int_pt, 0, 1, nn) = shp_mat;
+ }
+
+ // 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());
+
+ // TODO
+ // optimization: Store decomposition of N*N^T or N^T for reuse?
+ // cf. http://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html
+ // cf. http://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html
+ // cf. http://eigen.tuxfamily.org/dox/classEigen_1_1LDLT.html
+ // cf. https://eigen.tuxfamily.org/dox/classEigen_1_1HouseholderQR.html
+ // Extrapolate several values at once?
+
+ // do the least squares computation using QR decomposition.
+ Eigen::VectorXd tmp = N.householderQr().solve(integration_point_values_vec);
+
+ // option: do the least squares computation using LDLT decomposition.
+ // Eigen::VectorXd tmp =
+ // (N*N.transpose()).ldlt().solve(N*integration_point_values_vec);
+
+ // 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?
+ counts.add(global_indices, std::vector<double>(global_indices.size(), 1.0));
+}
+
+void LocalLinearLeastSquaresExtrapolator::calculateResidualElement(
+ std::size_t const element_index,
+ ExtrapolatableElementCollection const& extrapolatables)
+{
+ 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;
+
+ // 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) {
+ // TODO PETSc negative indices?
+ nodal_vals_element[i] = _nodal_values[global_indices[i]];
+ }
+
+ double residual = 0.0;
+ double gp_val_extrapol = 0.0;
+
+ for (unsigned gp = 0; gp < ni; ++gp) {
+ NumLib::shapeFunctionInterpolate(
+ 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;
+ }
+
+ _residuals.set(element_index, std::sqrt(residual / ni));
+}
+
+} // namespace NumLib
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index df3f0ffb093eb3ec788929846ff4e95cf72f7c85..52766bd8bda6b5c761a336876302d71845ab56ec 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
@@ -35,50 +35,20 @@ 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
+ /*! Constructs a new instance.
*
* \note
- * The \c dof_table of \c matrix_specs must be set, and it must point to a
- * dof_table for one single component variable.
+ * The \c dof_table must point to a d.o.f. table for one single-component
+ * variable.
*/
explicit LocalLinearLeastSquaresExtrapolator(
- MathLib::MatrixSpecifications const& matrix_specs,
- NumLib::LocalToGlobalIndexMap const& dof_table)
- : _nodal_values(
- NumLib::GlobalVectorProvider::provider.getVector(
- matrix_specs))
-#ifndef USE_PETSC
- ,
- _residuals(dof_table.size())
-#else
- ,
- _residuals(dof_table.size(), false)
-#endif
- ,
- _local_to_global(dof_table)
- {
- /* Note in case the following assertion fails.
- * If you copied the extrapolation code, for your processes from
- * somewhere, note that the code from the groundwater flow process might
- * not suit your needs: It is a special case and is therefore most
- * likely too simplistic. You better adapt the extrapolation code from
- * some more advanced process, like the TES process.
- */
- assert(dof_table.getNumberOfComponents() == 1 &&
- "The d.o.f. table passed must be for one variable that has "
- "only one component!");
- }
+ NumLib::LocalToGlobalIndexMap const& dof_table);
- void extrapolate(LocalAssemblers const& local_assemblers,
- PropertyTag const property) override;
+ void extrapolate(
+ ExtrapolatableElementCollection const& extrapolatables) override;
/*! \copydoc Extrapolator::calculateResiduals()
*
@@ -87,8 +57,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 +78,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 calculateResidualElement(
+ std::size_t const element_index,
+ ExtrapolatableElementCollection const& extrapolatables);
GlobalVector& _nodal_values; //!< extrapolated nodal values
GlobalVector _residuals; //!< extrapolation residuals
@@ -129,8 +100,7 @@ private:
//! Avoids frequent reallocations.
std::vector<double> _integration_point_values_cache;
};
-}
-#include "LocalLinearLeastSquaresExtrapolator-impl.h"
+} // namespace NumLib
#endif // NUMLIB_LOCAL_LLSQ_EXTRAPOLATOR_H
diff --git a/ProcessLib/ExtrapolatorData.h b/ProcessLib/ExtrapolatorData.h
new file mode 100644
index 0000000000000000000000000000000000000000..15e550b73ba423b049543b8ba8cbeb21d6c3ca4a
--- /dev/null
+++ b/ProcessLib/ExtrapolatorData.h
@@ -0,0 +1,102 @@
+/**
+ * \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 PROCESSLIB_EXTRAPOLATORDATA_H
+#define PROCESSLIB_EXTRAPOLATORDATA_H
+
+#include <memory>
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Extrapolation/Extrapolator.h"
+
+namespace ProcessLib
+{
+/*! Helper struct containing an extrapolator and a single component DOF table.
+ *
+ * Storage for the DOF table is managed optionally.
+ *
+ * \todo Later on this struct shall be moved, e.g., be merged with the process
+ * output class.
+ */
+class ExtrapolatorData
+{
+public:
+ ExtrapolatorData() = default;
+
+ /*! Constructs a new instance.
+ *
+ * \param extrapolator the extrapolator managed by the instance being
+ * created
+ * \param dof_table_single_component the d.o.f. table used by the \c
+ * extrapolator
+ * \param manage_storage If true the memory of \c dof_table_single_component
+ * will be freed by the destructor of this class.
+ */
+ ExtrapolatorData(
+ std::unique_ptr<NumLib::Extrapolator>&& extrapolator,
+ NumLib::LocalToGlobalIndexMap const* const dof_table_single_component,
+ bool const manage_storage)
+ : _extrapolator(std::move(extrapolator))
+ , _dof_table_single_component(dof_table_single_component)
+ , _manage_storage(manage_storage)
+ {
+ }
+
+ ExtrapolatorData(ExtrapolatorData&& other)
+ : _extrapolator(std::move(other._extrapolator))
+ , _dof_table_single_component(other._dof_table_single_component)
+ , _manage_storage(other._manage_storage)
+ {
+ other._manage_storage = false;
+ other._dof_table_single_component = nullptr;
+ }
+
+ ExtrapolatorData& operator=(ExtrapolatorData&& other)
+ {
+ cleanup();
+ _manage_storage = other._manage_storage;
+ _dof_table_single_component = other._dof_table_single_component;
+ _extrapolator = std::move(other._extrapolator);
+ other._dof_table_single_component = nullptr;
+ other._manage_storage = false;
+ return *this;
+ }
+
+ NumLib::LocalToGlobalIndexMap const& getDOFTable() const
+ {
+ return *_dof_table_single_component;
+ }
+ NumLib::Extrapolator& getExtrapolator() const { return *_extrapolator; }
+
+ ~ExtrapolatorData() { cleanup(); }
+
+private:
+ //! Deletes the d.o.f table if it is allowed to do so.
+ void cleanup()
+ {
+ if (_manage_storage) {
+ delete _dof_table_single_component;
+ _dof_table_single_component = nullptr;
+ }
+ }
+
+ //! Extrapolator managed by the ExtrapolatorData instance.
+ std::unique_ptr<NumLib::Extrapolator> _extrapolator;
+
+ //! D.o.f. table used by the extrapolator.
+ NumLib::LocalToGlobalIndexMap const* _dof_table_single_component = nullptr;
+
+ //! If true, free storage of the d.o.f. table in the ExtrapolatorData
+ //! destructor.
+ bool _manage_storage = false;
+};
+
+} // namespace ProcessLib
+
+
+#endif // PROCESSLIB_EXTRAPOLATORDATA_H
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 6f309631377062dae369d2b72ca48164944cbf04..0eaef58f200ba5215d6e0c50b57642600b06a262 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -12,6 +12,7 @@
#include <vector>
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "ProcessLib/LocalAssemblerInterface.h"
@@ -26,18 +27,22 @@ namespace ProcessLib
namespace GroundwaterFlow
{
-enum class IntegrationPointValue {
- DarcyVelocityX,
- DarcyVelocityY,
- DarcyVelocityZ
-};
-
const unsigned NUM_NODAL_DOF = 1;
class GroundwaterFlowLocalAssemblerInterface
: public ProcessLib::LocalAssemblerInterface
- , public NumLib::Extrapolatable<IntegrationPointValue>
-{};
+ , public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const = 0;
+};
template <typename ShapeFunction,
typename IntegrationMethod,
@@ -118,22 +123,24 @@ public:
}
std::vector<double> const&
- getIntegrationPointValues(IntegrationPointValue const property,
- std::vector<double>& /*cache*/) const override
+ getIntPtDarcyVelocityX(std::vector<double>& /*cache*/) const override
{
- switch (property)
- {
- case IntegrationPointValue::DarcyVelocityX:
- return _darcy_velocities[0];
- case IntegrationPointValue::DarcyVelocityY:
- assert(GlobalDim > 1);
- return _darcy_velocities[1];
- case IntegrationPointValue::DarcyVelocityZ:
- assert(GlobalDim > 2);
- return _darcy_velocities[2];
- }
+ assert(_darcy_velocities.size() > 1);
+ return _darcy_velocities[0];
+ }
+
+ std::vector<double> const&
+ getIntPtDarcyVelocityY(std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 1);
+ return _darcy_velocities[1];
+ }
- OGS_FATAL("");
+ std::vector<double> const&
+ getIntPtDarcyVelocityZ(std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 2);
+ return _darcy_velocities[2];
}
private:
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index e25c277069c31cd3c1e787efa2944df7be25888f..6aa0b2a37c7814eb3cb8f7468c6c62dec7d11a13 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -53,28 +53,25 @@ void GroundwaterFlowProcess::initializeConcreteProcess(
mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
_local_assemblers, _process_data);
- // TODO Later on the DOF table can change during the simulation!
- _extrapolator.reset(new ExtrapolatorImplementation(
- Base::getMatrixSpecifications(), *Base::_local_to_global_index_map));
-
- Base::_secondary_variables.addSecondaryVariable(
+ _secondary_variables.addSecondaryVariable(
"darcy_velocity_x", 1,
- makeExtrapolator(IntegrationPointValue::DarcyVelocityX, *_extrapolator,
- _local_assemblers));
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &GroundwaterFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
- if (mesh.getDimension() > 1)
- {
- Base::_secondary_variables.addSecondaryVariable(
+ if (mesh.getDimension() > 1) {
+ _secondary_variables.addSecondaryVariable(
"darcy_velocity_y", 1,
- makeExtrapolator(IntegrationPointValue::DarcyVelocityY,
- *_extrapolator, _local_assemblers));
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &GroundwaterFlowLocalAssemblerInterface::
+ getIntPtDarcyVelocityY));
}
- if (mesh.getDimension() > 2)
- {
- Base::_secondary_variables.addSecondaryVariable(
+ if (mesh.getDimension() > 2) {
+ _secondary_variables.addSecondaryVariable(
"darcy_velocity_z", 1,
- makeExtrapolator(IntegrationPointValue::DarcyVelocityZ,
- *_extrapolator, _local_assemblers));
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &GroundwaterFlowLocalAssemblerInterface::
+ getIntPtDarcyVelocityZ));
}
}
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index adaaf66d597dcaac9aee5336b5b892d4e445ee45..1c5254d2f2b708bf8af7c8f996155bc2bf22a8c5 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -42,13 +42,6 @@ public:
//! @}
private:
- using ExtrapolatorInterface =
- NumLib::Extrapolator<IntegrationPointValue,
- GroundwaterFlowLocalAssemblerInterface>;
- using ExtrapolatorImplementation =
- NumLib::LocalLinearLeastSquaresExtrapolator<
- IntegrationPointValue, GroundwaterFlowLocalAssemblerInterface>;
-
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh,
@@ -62,8 +55,6 @@ private:
std::vector<std::unique_ptr<GroundwaterFlowLocalAssemblerInterface>>
_local_assemblers;
-
- std::unique_ptr<ExtrapolatorInterface> _extrapolator;
};
} // namespace GroundwaterFlow
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 10258da6e888cbb9a8c5a9834583c641abbe6b03..f0cb228849450fca0d3041b3135d42ba7a644477 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -10,6 +10,7 @@
#include "Process.h"
#include "NumLib/DOF/ComputeSparsityPattern.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "DirichletBc.h"
#include "NeumannBc.h"
#include "NeumannBcAssembler.h"
@@ -52,6 +53,9 @@ void Process::initialize()
DBUG("Compute sparsity pattern");
computeSparsityPattern();
+ DBUG("Initialize the extrapolator");
+ initializeExtrapolator();
+
initializeConcreteProcess(*_local_to_global_index_map, _mesh,
_integration_order);
@@ -168,6 +172,42 @@ void Process::constructDofTable()
std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_LOCATION));
}
+void Process::initializeExtrapolator()
+{
+ NumLib::LocalToGlobalIndexMap const* dof_table_single_component;
+ bool manage_storage;
+
+ if (_local_to_global_index_map->getNumberOfComponents() == 1)
+ {
+ // For single-variable-single-component processes reuse the existing DOF
+ // table.
+ dof_table_single_component = _local_to_global_index_map.get();
+ manage_storage = false;
+ }
+ else
+ {
+ // Otherwise construct a new DOF table.
+ std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
+ all_mesh_subsets_single_component;
+ all_mesh_subsets_single_component.emplace_back(
+ new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get()));
+
+ dof_table_single_component = new NumLib::LocalToGlobalIndexMap(
+ std::move(all_mesh_subsets_single_component),
+ // by location order is needed for output
+ NumLib::ComponentOrder::BY_LOCATION);
+ manage_storage = true;
+ }
+
+ std::unique_ptr<NumLib::Extrapolator> extrapolator(
+ new NumLib::LocalLinearLeastSquaresExtrapolator(
+ *dof_table_single_component));
+
+ // TODO Later on the DOF table can change during the simulation!
+ _extrapolator_data = ExtrapolatorData(
+ std::move(extrapolator), dof_table_single_component, manage_storage);
+}
+
void Process::setInitialConditions(ProcessVariable const& variable,
int const variable_id,
int const component_id,
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 72d4009202b733b08518e5116292c273fcb5d3f9..ccedc3ae2df699f732db2e0d0e8e40d9c3fde951 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -14,6 +14,7 @@
#include "NumLib/ODESolver/ODESystem.h"
#include "NumLib/ODESolver/TimeDiscretization.h"
+#include "ExtrapolatorData.h"
#include "Parameter.h"
#include "ProcessOutput.h"
#include "SecondaryVariable.h"
@@ -85,6 +86,17 @@ public:
return *_time_discretization;
}
+protected:
+ NumLib::Extrapolator& getExtrapolator() const
+ {
+ return _extrapolator_data.getExtrapolator();
+ }
+
+ NumLib::LocalToGlobalIndexMap const& getSingleComponentDOFTable() const
+ {
+ return _extrapolator_data.getDOFTable();
+ }
+
private:
/// Process specific initialization called by initialize().
virtual void initializeConcreteProcess(
@@ -104,6 +116,8 @@ private:
void constructDofTable();
+ void initializeExtrapolator();
+
/// Sets the initial condition values in the solution vector x for a given
/// process variable and component.
void setInitialConditions(ProcessVariable const& variable,
@@ -142,6 +156,8 @@ private:
/// Variables used by this process.
std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
+
+ ExtrapolatorData _extrapolator_data;
};
/// Find process variables in \c variables whose names match the settings under
diff --git a/ProcessLib/SecondaryVariable.h b/ProcessLib/SecondaryVariable.h
index bf99e2e4c3fb4dfef1072c3be335b9caac102063..f845c1c244c5ac7f9fc9e1f4d2abd3c7c003204b 100644
--- a/ProcessLib/SecondaryVariable.h
+++ b/ProcessLib/SecondaryVariable.h
@@ -13,8 +13,12 @@
#include "BaseLib/ConfigTree.h"
#include "BaseLib/uniqueInsert.h"
#include "NumLib/Extrapolation/Extrapolator.h"
+#include "NumLib/Extrapolation/ExtrapolatableElementCollection.h"
-namespace NumLib { class LocalToGlobalIndexMap; }
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}
namespace ProcessLib
{
@@ -156,44 +160,50 @@ private:
std::set<std::string> _all_secondary_variables;
};
-
-//! Creates an object that computes a secondary variable via extrapolation
-//! of integration point values.
-template<typename PropertyEnum, typename LocalAssembler>
-SecondaryVariableFunctions
-makeExtrapolator(PropertyEnum const property,
- NumLib::Extrapolator<PropertyEnum, LocalAssembler>&
- extrapolator,
- typename NumLib::Extrapolator<PropertyEnum,
- LocalAssembler>::LocalAssemblers const& local_assemblers)
+/*! Creates an object that computes a secondary variable via extrapolation of
+ * integration point values.
+ *
+ * \param extrapolator The extrapolator used for extrapolation.
+ * \param local_assemblers The collection of local assemblers whose integration
+ * point values will be extrapolated.
+ * \param integration_point_values_method The member function of the local
+ * assembler returning/computing the integration point values of the specific
+ * property being extrapolated.
+ */
+template <typename LocalAssemblerCollection>
+SecondaryVariableFunctions makeExtrapolator(
+ NumLib::Extrapolator& extrapolator,
+ LocalAssemblerCollection const& local_assemblers,
+ typename NumLib::ExtrapolatableLocalAssemblerCollection<
+ LocalAssemblerCollection>::IntegrationPointValuesMethod
+ integration_point_values_method)
{
- static_assert(std::is_base_of<
- NumLib::Extrapolatable<PropertyEnum>, LocalAssembler>::value,
- "The passed local assembler type (i.e. the local assembler interface) must"
- " derive from NumLib::Extrapolatable<>.");
-
- auto const eval_field = [property, &extrapolator, &local_assemblers](
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::unique_ptr<GlobalVector>& /*result_cache*/
- ) -> GlobalVector const&
- {
- extrapolator.extrapolate(local_assemblers, property);
+ auto const eval_field = [&extrapolator, &local_assemblers,
+ integration_point_values_method](
+ GlobalVector const& /*x*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::unique_ptr<GlobalVector> & /*result_cache*/
+ ) -> GlobalVector const& {
+ auto const extrapolatables = NumLib::makeExtrapolatable(
+ local_assemblers, integration_point_values_method);
+ extrapolator.extrapolate(extrapolatables);
return extrapolator.getNodalValues();
};
- auto const eval_residuals = [property, &extrapolator, &local_assemblers](
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
- std::unique_ptr<GlobalVector>& /*result_cache*/
- ) -> GlobalVector const&
- {
- extrapolator.calculateResiduals(local_assemblers, property);
+ auto const eval_residuals = [&extrapolator, &local_assemblers,
+ integration_point_values_method](
+ GlobalVector const& /*x*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::unique_ptr<GlobalVector> & /*result_cache*/
+ ) -> GlobalVector const& {
+ auto const extrapolatables = NumLib::makeExtrapolatable(
+ local_assemblers, integration_point_values_method);
+ extrapolator.calculateResiduals(extrapolatables);
return extrapolator.getElementResiduals();
};
- return { eval_field, eval_residuals };
+ return {eval_field, eval_residuals};
}
-} // namespace ProcessLib
+} // namespace ProcessLib
#endif // PROCESSLIB_SECONDARY_VARIABLE_H
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 331f67b63523159a21b2fcf18215b9d90e4061e7..4eae782c384713e7674833344b99cfab7446d6c2 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -194,10 +194,81 @@ template <typename ShapeFunction_, typename IntegrationMethod_,
unsigned GlobalDim>
std::vector<double> const& TESLocalAssembler<
ShapeFunction_, IntegrationMethod_,
- GlobalDim>::getIntegrationPointValues(TESIntPtVariables const var,
- std::vector<double>& cache) const
+ GlobalDim>::getIntPtSolidDensity(std::vector<double>& /*cache*/) const
{
- return _d.getIntegrationPointValues(var, cache);
+ return _d.getData().solid_density;
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_,
+ GlobalDim>::getIntPtLoading(std::vector<double>& cache) const
+{
+ auto const rho_SR = _d.getData().solid_density;
+ auto const rho_SR_dry = _d.getAssemblyParameters().rho_SR_dry;
+
+ cache.clear();
+ cache.reserve(rho_SR.size());
+
+ for (auto const rho : rho_SR) {
+ cache.push_back(rho/rho_SR_dry - 1.0);
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ unsigned GlobalDim>
+std::vector<double> const&
+TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+ getIntPtReactionDampingFactor(std::vector<double>& cache) const
+{
+ auto const fac = _d.getData().reaction_adaptor->getReactionDampingFactor();
+ auto const num_integration_points = _d.getData().solid_density.size();
+
+ cache.clear();
+ cache.resize(num_integration_points, fac);
+
+ return cache;
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_,
+ GlobalDim>::getIntPtReactionRate(std::vector<double>& /*cache*/) const
+{
+ return _d.getData().reaction_rate;
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_,
+ GlobalDim>::getIntPtDarcyVelocityX(std::vector<double>& /*cache*/) const
+{
+ return _d.getData().velocity[0];
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_,
+ GlobalDim>::getIntPtDarcyVelocityY(std::vector<double>& /*cache*/) const
+{
+ assert(_d.getData().velocity.size() > 1);
+ return _d.getData().velocity[1];
+}
+
+template <typename ShapeFunction_, typename IntegrationMethod_,
+ unsigned GlobalDim>
+std::vector<double> const& TESLocalAssembler<
+ ShapeFunction_, IntegrationMethod_,
+ GlobalDim>::getIntPtDarcyVelocityZ(std::vector<double>& /*cache*/) const
+{
+ assert(_d.getData().velocity.size() > 2);
+ return _d.getData().velocity[2];
}
template <typename ShapeFunction_, typename IntegrationMethod_,
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index 37ce25c30459fb2f945a8eac62ec0086c0589fb6..1826b2c6f4cc2fec5448b5559c80aa37d3c914b5 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -14,24 +14,44 @@
#include <vector>
#include "ProcessLib/LocalAssemblerInterface.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "TESAssemblyParams.h"
#include "TESLocalAssemblerInner-fwd.h"
-#include "NumLib/Extrapolation/Extrapolator.h"
-
namespace ProcessLib
{
namespace TES
{
class TESLocalAssemblerInterface
: public ProcessLib::LocalAssemblerInterface,
- public NumLib::Extrapolatable<TESIntPtVariables>
+ public NumLib::ExtrapolatableElement
{
public:
virtual ~TESLocalAssemblerInterface() = default;
virtual bool checkBounds(std::vector<double> const& local_x,
std::vector<double> const& local_x_prev_ts) = 0;
+
+ virtual std::vector<double> const& getIntPtSolidDensity(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtLoading(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtReactionDampingFactor(
+ std::vector<double>& cache) const = 0;
+
+ virtual std::vector<double> const& getIntPtReactionRate(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const = 0;
};
template <typename ShapeFunction_, typename IntegrationMethod_,
@@ -66,9 +86,26 @@ public:
bool checkBounds(std::vector<double> const& local_x,
std::vector<double> const& local_x_prev_ts) override;
- std::vector<double> const& getIntegrationPointValues(
- TESIntPtVariables const var, std::vector<double>& cache) const override;
+ std::vector<double> const& getIntPtSolidDensity(
+ std::vector<double>& /*cache*/) const override;
+
+ std::vector<double> const& getIntPtLoading(
+ std::vector<double>& cache) const override;
+
+ std::vector<double> const& getIntPtReactionDampingFactor(
+ std::vector<double>& cache) const override;
+
+ std::vector<double> const& getIntPtReactionRate(
+ std::vector<double>& /*cache*/) const override;
+
+ std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const override;
+
+ std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const override;
+ std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const override;
private:
std::vector<ShapeMatrices> _shape_matrices;
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl.h b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
index d587920df981851367fd1f4e838c8fb2412595d5..b0d63166cd60ad493793a48d369cc3de3a234d52 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner-impl.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
@@ -218,57 +218,6 @@ void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
_d.rho_GR = fluid_density(_d.p, _d.T, _d.vapour_mass_fraction);
}
-template <typename Traits>
-std::vector<double> const&
-TESLocalAssemblerInner<Traits>::getIntegrationPointValues(
- TESIntPtVariables const var, std::vector<double>& cache) const
-{
- switch (var)
- {
- case TESIntPtVariables::REACTION_RATE:
- return _d.reaction_rate;
- case TESIntPtVariables::SOLID_DENSITY:
- return _d.solid_density;
- case TESIntPtVariables::VELOCITY_X:
- return _d.velocity[0];
- case TESIntPtVariables::VELOCITY_Y:
- assert(_d.velocity.size() >= 2);
- return _d.velocity[1];
- case TESIntPtVariables::VELOCITY_Z:
- assert(_d.velocity.size() >= 3);
- return _d.velocity[2];
-
- case TESIntPtVariables::LOADING:
- {
- auto& Cs = cache;
- Cs.clear();
- Cs.reserve(_d.solid_density.size());
-
- for (auto rho_SR : _d.solid_density)
- {
- Cs.push_back(rho_SR / _d.ap.rho_SR_dry - 1.0);
- }
-
- return Cs;
- }
-
- // TODO that's an element value, ain't it?
- case TESIntPtVariables::REACTION_DAMPING_FACTOR:
- {
- auto const num_integration_points = _d.solid_density.size();
- auto& alphas = cache;
- alphas.clear();
- alphas.resize(num_integration_points,
- _d.reaction_adaptor->getReactionDampingFactor());
-
- return alphas;
- }
- }
-
- cache.clear();
- return cache;
-}
-
template <typename Traits>
void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
unsigned integration_point,
diff --git a/ProcessLib/TES/TESLocalAssemblerInner.h b/ProcessLib/TES/TESLocalAssemblerInner.h
index 3e36cb4718059fb80bb4ffd989008ced34fcb05d..77d4638fc507c22d697a7cc89dee6f0bc5b39587 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner.h
@@ -23,17 +23,6 @@ namespace ProcessLib
{
namespace TES
{
-enum class TESIntPtVariables : unsigned
-{
- SOLID_DENSITY,
- REACTION_RATE,
- VELOCITY_X,
- VELOCITY_Y,
- VELOCITY_Z,
- LOADING,
- REACTION_DAMPING_FACTOR
-};
-
template <typename Traits>
class TESLocalAssemblerInner
{
@@ -56,9 +45,6 @@ public:
void preEachAssemble();
- std::vector<double> const& getIntegrationPointValues(
- TESIntPtVariables const var, std::vector<double>& cache) const;
-
// TODO better encapsulation
AssemblyParams const& getAssemblyParameters() const { return _d.ap; }
TESFEMReactionAdaptor const& getReactionAdaptor() const
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index b35eb3c189645aa58fc0ca75b8874d2bc1bbbd69..f2c373a56acf549d040aabda67c71392b6906d6b 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -175,28 +175,6 @@ void TESProcess::initializeConcreteProcess(
mesh.getDimension(), mesh.getElements(), dof_table, integration_order,
_local_assemblers, _assembly_params);
- // TODO move the two data members somewhere else.
- // for extrapolation of secondary variables
- std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
- all_mesh_subsets_single_component;
- all_mesh_subsets_single_component.emplace_back(
- new MeshLib::MeshSubsets(this->_mesh_subset_all_nodes.get()));
- _local_to_global_index_map_single_component.reset(
- new NumLib::LocalToGlobalIndexMap(
- std::move(all_mesh_subsets_single_component),
- // by location order is needed for output
- NumLib::ComponentOrder::BY_LOCATION));
-
- {
- auto const& l = *_local_to_global_index_map_single_component;
- _extrapolator.reset(new ExtrapolatorImplementation(
- MathLib::MatrixSpecifications(l.dofSizeWithoutGhosts(),
- l.dofSizeWithoutGhosts(),
- &l.getGhostIndices(),
- nullptr),
- l));
- }
-
// secondary variables
auto add2nd = [&](std::string const& var_name, unsigned const n_components,
SecondaryVariableFunctions&& fcts) {
@@ -204,22 +182,29 @@ void TESProcess::initializeConcreteProcess(
std::move(fcts));
};
auto makeEx =
- [&](TESIntPtVariables var) -> SecondaryVariableFunctions {
- return ProcessLib::makeExtrapolator(var, *_extrapolator,
- _local_assemblers);
+ [&](std::vector<double> const& (TESLocalAssemblerInterface::*method)(
+ std::vector<double>&)const) -> SecondaryVariableFunctions {
+ return ProcessLib::makeExtrapolator(getExtrapolator(),
+ _local_assemblers, method);
};
- add2nd("solid_density", 1, makeEx(TESIntPtVariables::SOLID_DENSITY));
- add2nd("reaction_rate", 1, makeEx(TESIntPtVariables::REACTION_RATE));
- add2nd("velocity_x", 1, makeEx(TESIntPtVariables::VELOCITY_X));
+ add2nd("solid_density", 1,
+ makeEx(&TESLocalAssemblerInterface::getIntPtSolidDensity));
+ add2nd("reaction_rate", 1,
+ makeEx(&TESLocalAssemblerInterface::getIntPtReactionRate));
+
+ add2nd("velocity_x", 1,
+ makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityX));
if (mesh.getDimension() >= 2)
- add2nd("velocity_y", 1, makeEx(TESIntPtVariables::VELOCITY_Y));
+ add2nd("velocity_y", 1,
+ makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityY));
if (mesh.getDimension() >= 3)
- add2nd("velocity_z", 1, makeEx(TESIntPtVariables::VELOCITY_Z));
+ add2nd("velocity_z", 1,
+ makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityZ));
- add2nd("loading", 1, makeEx(TESIntPtVariables::LOADING));
+ add2nd("loading", 1, makeEx(&TESLocalAssemblerInterface::getIntPtLoading));
add2nd("reaction_damping_factor", 1,
- makeEx(TESIntPtVariables::REACTION_DAMPING_FACTOR));
+ makeEx(&TESLocalAssemblerInterface::getIntPtReactionDampingFactor));
namespace PH = std::placeholders;
using Self = TESProcess;
@@ -338,7 +323,7 @@ TESProcess::computeVapourPartialPressure(
{
assert(&dof_table == this->_local_to_global_index_map.get());
- auto const& dof_table_single = *_local_to_global_index_map_single_component;
+ auto const& dof_table_single = getSingleComponentDOFTable();
result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
{dof_table_single.dofSizeWithoutGhosts(),
dof_table_single.dofSizeWithoutGhosts(),
@@ -371,7 +356,7 @@ TESProcess::computeRelativeHumidity(
{
assert(&dof_table == this->_local_to_global_index_map.get());
- auto const& dof_table_single = *_local_to_global_index_map_single_component;
+ auto const& dof_table_single = getSingleComponentDOFTable();
result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
{dof_table_single.dofSizeWithoutGhosts(),
dof_table_single.dofSizeWithoutGhosts(),
@@ -409,7 +394,7 @@ TESProcess::computeEquilibriumLoading(
{
assert(&dof_table == this->_local_to_global_index_map.get());
- auto const& dof_table_single = *_local_to_global_index_map_single_component;
+ auto const& dof_table_single = getSingleComponentDOFTable();
result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
{dof_table_single.dofSizeWithoutGhosts(),
dof_table_single.dofSizeWithoutGhosts(),
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 987f84be056ade216626e93881447609a1217989..ec60ca4c03c5441557253149e5eb70a65037216a 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -31,8 +31,6 @@ namespace TES
{
class TESProcess final : public Process
{
- using BP = Process; //!< "Base Process"
-
public:
TESProcess(
MeshLib::Mesh& mesh,
@@ -51,13 +49,8 @@ public:
NumLib::IterationResult postIteration(GlobalVector const& x) override;
bool isLinear() const override { return false; }
-private:
- using ExtrapolatorInterface =
- NumLib::Extrapolator<TESIntPtVariables, TESLocalAssemblerInterface>;
- using ExtrapolatorImplementation =
- NumLib::LocalLinearLeastSquaresExtrapolator<
- TESIntPtVariables, TESLocalAssemblerInterface>;
+private:
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh, unsigned const integration_order) override;
@@ -85,11 +78,6 @@ private:
AssemblyParams _assembly_params;
- std::unique_ptr<NumLib::LocalToGlobalIndexMap>
- _local_to_global_index_map_single_component;
-
- std::unique_ptr<ExtrapolatorInterface> _extrapolator;
-
// used for checkBounds()
std::unique_ptr<GlobalVector> _x_previous_timestep;
};
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index b76cb3b76847dccf38d8479c4feab4298f2ee0f4..c33916bfc9c807daa1c6d9d33c388c1ce241e277 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -10,26 +10,26 @@
#include <random>
#include <gtest/gtest.h>
-#include "NumLib/DOF/MatrixProviderUser.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
-
#include "MathLib/LinAlg/LinAlg.h"
#include "MeshLib/MeshGenerators/MeshGenerator.h"
#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/DOF/MatrixProviderUser.h"
+#include "NumLib/Extrapolation/ExtrapolatableElementCollection.h"
#include "NumLib/Extrapolation/Extrapolator.h"
#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
-
+#include "NumLib/Function/Interpolation.h"
#include "NumLib/NumericsConfig.h"
+
#include "ProcessLib/Utils/LocalDataInitializer.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
#include "ProcessLib/Utils/InitShapeMatrices.h"
-
namespace
{
@@ -63,26 +63,25 @@ void fillVectorRandomly(Vector& x)
}
-enum class IntegrationPointValue
-{
- StoredQuantity, // some quantity acutally stored in the local assembler
- DerivedQuantity // a quantity computed for each integration point on-the-fly
-};
-
-class LocalAssemblerDataInterface
- : public NumLib::Extrapolatable<IntegrationPointValue>
+class LocalAssemblerDataInterface : public NumLib::ExtrapolatableElement
{
public:
virtual void interpolateNodalValuesToIntegrationPoints(
- std::vector<double> const& local_nodal_values,
- IntegrationPointValue const property) = 0;
+ std::vector<double> const& local_nodal_values) = 0;
+
+ virtual std::vector<double> const& getStoredQuantity(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getDerivedQuantity(
+ std::vector<double>& cache) const = 0;
};
-template<typename ShapeFunction,
- typename IntegrationMethod,
- unsigned GlobalDim>
-class LocalAssemblerData
- : public LocalAssemblerDataInterface
+using IntegrationPointValuesMethod = std::vector<double> const& (
+ LocalAssemblerDataInterface::*)(std::vector<double>&)const;
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class LocalAssemblerData : public LocalAssemblerDataInterface
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
@@ -93,9 +92,11 @@ public:
unsigned const integration_order)
: _shape_matrices(
ProcessLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
- IntegrationMethod, GlobalDim>(e, integration_order))
+ IntegrationMethod, GlobalDim>(
+ e, integration_order))
, _int_pt_values(_shape_matrices.size())
- {}
+ {
+ }
Eigen::Map<const Eigen::RowVectorXd>
getShapeMatrix(const unsigned integration_point) const override
@@ -106,42 +107,31 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
- std::vector<double> const&
- getIntegrationPointValues(IntegrationPointValue const property,
- std::vector<double>& cache) const override
+ std::vector<double> const& getStoredQuantity(
+ std::vector<double>& /*cache*/) const override
{
- switch (property)
- {
- case IntegrationPointValue::StoredQuantity:
- return _int_pt_values;
- case IntegrationPointValue::DerivedQuantity:
- cache.clear();
- for (auto value : _int_pt_values)
- cache.push_back(2.0 * value);
- return cache;
- }
+ return _int_pt_values;
+ }
- OGS_FATAL("");
+ std::vector<double> const& getDerivedQuantity(
+ std::vector<double>& cache) const override
+ {
+ cache.clear();
+ for (auto value : _int_pt_values)
+ cache.push_back(2.0 * value);
+ return cache;
}
void interpolateNodalValuesToIntegrationPoints(
- std::vector<double> const& local_nodal_values,
- IntegrationPointValue const property) override
+ std::vector<double> const& local_nodal_values) override
{
- switch (property)
- {
- case IntegrationPointValue::StoredQuantity:
- ::interpolateNodalValuesToIntegrationPoints(
- local_nodal_values, _shape_matrices, _int_pt_values);
- break;
- case IntegrationPointValue::DerivedQuantity:
- break;
- }
+ ::interpolateNodalValuesToIntegrationPoints(
+ local_nodal_values, _shape_matrices, _int_pt_values);
}
private:
std::vector<ShapeMatrices> _shape_matrices;
- std::vector<double> _int_pt_values;
+ std::vector<double> _int_pt_values;
};
class TestProcess
@@ -149,11 +139,9 @@ class TestProcess
public:
using LocalAssembler = LocalAssemblerDataInterface;
- using ExtrapolatorInterface =
- NumLib::Extrapolator<IntegrationPointValue, LocalAssembler>;
+ using ExtrapolatorInterface = NumLib::Extrapolator;
using ExtrapolatorImplementation =
- NumLib::LocalLinearLeastSquaresExtrapolator<
- IntegrationPointValue, LocalAssembler>;
+ NumLib::LocalLinearLeastSquaresExtrapolator;
TestProcess(MeshLib::Mesh const& mesh, unsigned const integration_order)
: _integration_order(integration_order)
@@ -161,51 +149,44 @@ public:
{
std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
all_mesh_subsets.emplace_back(
- new MeshLib::MeshSubsets{&_mesh_subset_all_nodes});
+ new MeshLib::MeshSubsets{&_mesh_subset_all_nodes});
_dof_table.reset(new NumLib::LocalToGlobalIndexMap(
- std::move(all_mesh_subsets),
- NumLib::ComponentOrder::BY_COMPONENT));
+ std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_COMPONENT));
// Passing _dof_table works, because this process has only one variable
// and the variable has exactly one component.
- _extrapolator.reset(new ExtrapolatorImplementation(
- MathLib::MatrixSpecifications{_dof_table->dofSizeWithoutGhosts(),
- _dof_table->dofSizeWithoutGhosts(),
- &_dof_table->getGhostIndices(),
- nullptr},
- *_dof_table));
+ _extrapolator.reset(new ExtrapolatorImplementation(*_dof_table));
createAssemblers(mesh);
}
void interpolateNodalValuesToIntegrationPoints(
- GlobalVector const& global_nodal_values,
- IntegrationPointValue const property) const
+ GlobalVector const& global_nodal_values) const
{
auto cb = [](std::size_t id, LocalAssembler& loc_asm,
NumLib::LocalToGlobalIndexMap const& dof_table,
- GlobalVector const& x,
- IntegrationPointValue const property) {
+ GlobalVector const& x) {
auto const indices = NumLib::getIndices(id, dof_table);
auto const local_x = x.get(indices);
- loc_asm.interpolateNodalValuesToIntegrationPoints(local_x,
- property);
+ loc_asm.interpolateNodalValuesToIntegrationPoints(local_x);
};
GlobalExecutor::executeDereferenced(cb, _local_assemblers, *_dof_table,
- global_nodal_values, property);
+ global_nodal_values);
}
- std::pair<GlobalVector const*, GlobalVector const*>
- extrapolate(IntegrationPointValue const property) const
+ std::pair<GlobalVector const*, GlobalVector const*> extrapolate(
+ IntegrationPointValuesMethod method) const
{
- _extrapolator->extrapolate(_local_assemblers, property);
- _extrapolator->calculateResiduals(_local_assemblers, property);
+ auto const extrapolatables =
+ NumLib::makeExtrapolatable(_local_assemblers, method);
+ _extrapolator->extrapolate(extrapolatables);
+ _extrapolator->calculateResiduals(extrapolatables);
- return { &_extrapolator->getNodalValues(),
- &_extrapolator->getElementResiduals() };
+ return {&_extrapolator->getNodalValues(),
+ &_extrapolator->getElementResiduals()};
}
private:
@@ -249,11 +230,8 @@ private:
std::unique_ptr<ExtrapolatorInterface> _extrapolator;
};
-
-void extrapolate(TestProcess const& pcs,
- IntegrationPointValue property,
- GlobalVector const&
- expected_extrapolated_global_nodal_values,
+void extrapolate(TestProcess const& pcs, IntegrationPointValuesMethod method,
+ GlobalVector const& expected_extrapolated_global_nodal_values,
std::size_t const nnodes, std::size_t const nelements)
{
namespace LinAlg = MathLib::LinAlg;
@@ -261,8 +239,8 @@ void extrapolate(TestProcess const& pcs,
auto const tolerance_dx = 20.0 * std::numeric_limits<double>::epsilon();
auto const tolerance_res = 5.0 * std::numeric_limits<double>::epsilon();
- auto const result = pcs.extrapolate(property);
- auto const& x_extra = *result.first;
+ auto const result = pcs.extrapolate(method);
+ auto const& x_extra = *result.first;
auto const& residual = *result.second;
ASSERT_EQ(nnodes, x_extra.size());
@@ -322,20 +300,20 @@ TEST(NumLib, DISABLED_Extrapolation)
fillVectorRandomly(*x);
- pcs.interpolateNodalValuesToIntegrationPoints(
- *x, IntegrationPointValue::StoredQuantity);
+ pcs.interpolateNodalValuesToIntegrationPoints(*x);
- // test extrapolation of a quantity that is stored in the local assembler
- extrapolate(pcs, IntegrationPointValue::StoredQuantity,
- *x, nnodes, nelements);
+ // test extrapolation of a quantity that is stored in the local
+ // assembler
+ extrapolate(pcs, &LocalAssemblerDataInterface::getStoredQuantity, *x,
+ nnodes, nelements);
// expect 2*x as extraplation result for derived quantity
auto two_x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(*x);
LinAlg::axpy(*two_x, 1.0, *x); // two_x = x + x
- // test extrapolation of a quantity that is derived from some integration
- // point values
- extrapolate(pcs, IntegrationPointValue::DerivedQuantity,
+ // test extrapolation of a quantity that is derived from some
+ // integration point values
+ extrapolate(pcs, &LocalAssemblerDataInterface::getDerivedQuantity,
*two_x, nnodes, nelements);
}
}