diff --git a/NumLib/Extrapolation/ExtrapolatableElementCollection.h b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
index c07d719a79a71a8cfd4e3c7d1bb53f0312c9f582..ab47613cad507e716182f93d212450cd5a85341c 100644
--- a/NumLib/Extrapolation/ExtrapolatableElementCollection.h
+++ b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
@@ -9,11 +9,17 @@
#pragma once
+#include <functional>
#include <vector>
+
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
+
#include "ExtrapolatableElement.h"
namespace NumLib
{
+class LocalToGlobalIndexMap;
+
/*! Adaptor to get information needed by an Extrapolator from an "arbitrary"
* collection of elements (e.g., local assemblers).
*
@@ -50,7 +56,10 @@ public:
* \endparblock
*/
virtual std::vector<double> const& getIntegrationPointValues(
- std::size_t const id, std::vector<double>& cache) const = 0;
+ std::size_t const id, const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const = 0;
//! Returns the number of elements whose properties shall be extrapolated.
virtual std::size_t size() const = 0;
@@ -64,9 +73,9 @@ class ExtrapolatableLocalAssemblerCollection
{
public:
//! LocalAssemblerCollection contains many LocalAssembler's.
- using LocalAssembler = typename std::decay<decltype(
- *std::declval<LocalAssemblerCollection>()[static_cast<std::size_t>(
- 0)])>::type;
+ 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 "
@@ -80,8 +89,12 @@ public:
* For further information about the \c cache parameter see
* ExtrapolatableElementCollection::getIntegrationPointValues().
*/
- using IntegrationPointValuesMethod = std::vector<double> const& (
- LocalAssembler::*)(std::vector<double>& cache) const;
+ using IntegrationPointValuesMethod =
+ std::function<std::vector<double> const&(
+ LocalAssembler const& loc_asm, const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache)>;
/*! Constructs a new instance.
*
@@ -92,9 +105,9 @@ public:
*/
ExtrapolatableLocalAssemblerCollection(
LocalAssemblerCollection const& local_assemblers,
- IntegrationPointValuesMethod integration_point_values_method)
- : _local_assemblers(local_assemblers)
- , _integration_point_values_method(integration_point_values_method)
+ IntegrationPointValuesMethod const& integration_point_values_method)
+ : _local_assemblers(local_assemblers),
+ _integration_point_values_method{integration_point_values_method}
{
}
@@ -106,10 +119,14 @@ public:
}
std::vector<double> const& getIntegrationPointValues(
- std::size_t const id, std::vector<double>& cache) const override
+ std::size_t const id, const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const override
{
auto const& loc_asm = *_local_assemblers[id];
- return (loc_asm.*_integration_point_values_method)(cache);
+ return _integration_point_values_method(loc_asm, t, current_solution,
+ dof_table, cache);
}
std::size_t size() const override { return _local_assemblers.size(); }
diff --git a/NumLib/Extrapolation/Extrapolator.h b/NumLib/Extrapolation/Extrapolator.h
index fa13cace48383ba255688c0780dc9673a7c84bf8..c658057c146b5e9edaac52d37edde3fe7153ea60 100644
--- a/NumLib/Extrapolation/Extrapolator.h
+++ b/NumLib/Extrapolation/Extrapolator.h
@@ -17,6 +17,7 @@
namespace NumLib
{
+class LocalToGlobalIndexMap;
//! Interface for classes that extrapolate integration point values to nodal
//! values.
@@ -25,7 +26,11 @@ class Extrapolator
public:
//! Extrapolates the given \c property from the given local assemblers.
virtual void extrapolate(
- ExtrapolatableElementCollection const& extrapolatables) = 0;
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table) = 0;
/*! Computes residuals from the extrapolation of the given \c property.
*
@@ -34,7 +39,11 @@ public:
* \pre extrapolate() must have been called before with the same arguments.
*/
virtual void calculateResiduals(
- ExtrapolatableElementCollection const& extrapolatables) = 0;
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table) = 0;
//! Returns the extrapolated nodal values.
//! \todo Maybe write directly to a MeshProperty.
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
index 29c973ae10a4dc8ba548f89b990a0bb07a63d44d..6a59ef02d0e0b108afba924986dfc83dc7249475 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
@@ -23,18 +23,7 @@ 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.dofSizeWithoutGhosts(),
- dof_table.getGhostIndices(), false),
-#endif
- _local_to_global(dof_table)
+ : _dof_table_single_component(dof_table)
{
/* Note in case the following assertion fails:
* If you copied the extrapolation code, for your processes from
@@ -43,58 +32,112 @@ LocalLinearLeastSquaresExtrapolator::LocalLinearLeastSquaresExtrapolator(
* 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!");
+ if (dof_table.getNumberOfComponents() != 1)
+ OGS_FATAL(
+ "The d.o.f. table passed must be for one variable that has "
+ "only one component!");
}
void LocalLinearLeastSquaresExtrapolator::extrapolate(
- ExtrapolatableElementCollection const& extrapolatables)
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table)
{
- _nodal_values.setZero();
+ auto const num_nodal_dof_result =
+ _dof_table_single_component.dofSizeWithoutGhosts() * num_components;
+ if (!_nodal_values ||
+ static_cast<std::size_t>(_nodal_values->size()) != num_nodal_dof_result)
+ {
+ _nodal_values = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+ MathLib::MatrixSpecifications{
+ num_nodal_dof_result, num_nodal_dof_result, nullptr, nullptr});
+ }
+ _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);
+ 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);
+ for (std::size_t i = 0; i < size; ++i)
+ {
+ extrapolateElement(i, num_components, extrapolatables, t,
+ current_solution, dof_table, *counts);
}
- MathLib::LinAlg::componentwiseDivide(_nodal_values, _nodal_values, *counts);
+ MathLib::LinAlg::componentwiseDivide(*_nodal_values, *_nodal_values,
+ *counts);
}
void LocalLinearLeastSquaresExtrapolator::calculateResiduals(
- ExtrapolatableElementCollection const& extrapolatables)
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table)
{
- assert(static_cast<std::size_t>(_residuals.size()) ==
- extrapolatables.size());
+ auto const num_element_dof_result = static_cast<GlobalIndexType>(
+ _dof_table_single_component.size() * num_components);
+
+ if (!_residuals || _residuals->size() != num_element_dof_result)
+ {
+#ifndef USE_PETSC
+ _residuals.reset(new GlobalVector{num_element_dof_result});
+#else
+ _residuals.reset(new GlobalVector{num_element_dof_result, false});
+#endif
+ }
+
+ if (static_cast<std::size_t>(num_element_dof_result) !=
+ extrapolatables.size() * num_components)
+ {
+ OGS_FATAL("mismatch in number of D.o.F.");
+ }
auto const size = extrapolatables.size();
- for (std::size_t i=0; i<size; ++i) {
- calculateResidualElement(i, extrapolatables);
+ for (std::size_t i = 0; i < size; ++i)
+ {
+ calculateResidualElement(i, num_components, extrapolatables, t,
+ current_solution, dof_table);
}
}
void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
std::size_t const element_index,
+ const unsigned num_components,
ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table,
GlobalVector& counts)
{
auto const& integration_point_values =
extrapolatables.getIntegrationPointValues(
- element_index, _integration_point_values_cache);
+ element_index, t, current_solution, dof_table,
+ _integration_point_values_cache);
auto const& N_0 = extrapolatables.getShapeMatrix(element_index, 0);
- const auto num_nodes = static_cast<unsigned>(N_0.cols());
- const unsigned num_int_pts = integration_point_values.size();
+ const unsigned num_nodes = N_0.cols();
+ const unsigned num_values = integration_point_values.size();
+
+ if (num_values % num_components != 0)
+ OGS_FATAL(
+ "The number of computed integration point values is not divisable "
+ "by the number of num_components. Maybe the computed property is "
+ "not a %d-component vector for each integration point.",
+ num_components);
- assert(num_int_pts >= num_nodes &&
- "Least squares is not possible if there are more nodes than"
- "integration points.");
+ // number of integration points in the element
+ const auto num_int_pts = num_values / num_components;
+
+ if (num_int_pts < static_cast<decltype(num_int_pts)>(num_nodes))
+ OGS_FATAL(
+ "Least squares is not possible if there are more nodes than"
+ "integration points.");
auto const pair_it_inserted = _qr_decomposition_cache.emplace(
std::make_pair(num_nodes, num_int_pts), CachedData{});
@@ -146,51 +189,108 @@ void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
OGS_FATAL("The cached and the passed shapematrices differ.");
}
- // TODO make gp_vals an Eigen::VectorXd const& ?
- auto const integration_point_values_vec =
- MathLib::toVector(integration_point_values);
-
- // Apply the pre-computed pseudo-inverse.
- Eigen::VectorXd const nodal_values =
- cached_data.A_pinv * integration_point_values_vec;
-
// TODO: for now always zeroth component is used. This has to be extended if
// multi-component properties shall be extrapolated
- auto const& global_indices = _local_to_global(element_index, 0).rows;
+ auto const& global_indices =
+ _dof_table_single_component(element_index, 0).rows;
+
+ if (num_components == 1)
+ {
+ // TODO make gp_vals an Eigen::VectorXd const& ?
+ auto const integration_point_values_vec =
+ MathLib::toVector(integration_point_values);
+
+ // Apply the pre-computed pseudo-inverse.
+ Eigen::VectorXd const nodal_values =
+ cached_data.A_pinv * integration_point_values_vec;
+
+ // TODO does that give rise to PETSc problems?
+ _nodal_values->add(global_indices, nodal_values);
+ counts.add(global_indices,
+ std::vector<double>(global_indices.size(), 1.0));
+ }
+ else
+ {
+ // TODO make gp_vals an Eigen::VectorXd const& ?
+ auto const integration_point_values_mat = MathLib::toMatrix(
+ integration_point_values, num_components, num_int_pts);
+
+ // Apply the pre-computed pseudo-inverse.
+ Eigen::MatrixXd const nodal_values =
+ cached_data.A_pinv * integration_point_values_mat.transpose();
+
+ std::vector<GlobalIndexType> indices;
+ indices.reserve(num_components * global_indices.size());
+
+ // component-wise in nodal_values_flat, location-wise ordering in
+ // _nodal_values
+ for (unsigned comp = 0; comp < num_components; ++comp)
+ {
+ for (auto i : global_indices)
+ {
+ indices.push_back(num_components * i + comp);
+ }
+ }
- // TODO does that give rise to PETSc problems?
- _nodal_values.add(global_indices, nodal_values);
- counts.add(global_indices, std::vector<double>(global_indices.size(), 1.0));
+ // Nodal_values are passed as a raw pointer, because PETScVector and
+ // EigenVector implementations differ slightly.
+ _nodal_values->add(indices, nodal_values.data());
+ counts.add(indices, std::vector<double>(indices.size(), 1.0));
+ }
}
void LocalLinearLeastSquaresExtrapolator::calculateResidualElement(
std::size_t const element_index,
- ExtrapolatableElementCollection const& extrapolatables)
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table)
{
auto const& int_pt_vals = extrapolatables.getIntegrationPointValues(
- element_index, _integration_point_values_cache);
+ element_index, t, current_solution, dof_table,
+ _integration_point_values_cache);
- // TODO: for now always zeroth component is used
- const auto& global_indices = _local_to_global(element_index, 0).rows;
+ const unsigned num_values = int_pt_vals.size();
+ if (num_values % num_components != 0)
+ OGS_FATAL(
+ "The number of computed integration point values is not divisable "
+ "by the number of num_components. Maybe the computed property is "
+ "not a %d-component vector for each integration point.",
+ num_components);
- const unsigned num_int_pts = int_pt_vals.size();
- const unsigned num_nodes = global_indices.size();
+ // number of integration points in the element
+ const auto num_int_pts = num_values / num_components;
- // filter nodal values of the current element
- Eigen::VectorXd nodal_vals_element(num_nodes);
- for (unsigned i = 0; i < num_nodes; ++i) {
- // TODO PETSc negative indices?
- nodal_vals_element[i] = _nodal_values[global_indices[i]];
- }
+ // TODO: for now always zeroth component is used
+ const auto& global_indices =
+ _dof_table_single_component(element_index, 0).rows;
+ const unsigned num_nodes = global_indices.size();
auto const& interpolation_matrix =
_qr_decomposition_cache.find({num_nodes, num_int_pts})->second.A;
- double const residual = (interpolation_matrix * nodal_vals_element -
- MathLib::toVector(int_pt_vals))
- .squaredNorm();
+ Eigen::VectorXd nodal_vals_element(num_nodes);
+ auto const int_pt_vals_mat =
+ MathLib::toMatrix(int_pt_vals, num_components, num_int_pts);
+
+ for (unsigned comp = 0; comp < num_components; ++comp)
+ {
+ // filter nodal values of the current element
+ for (unsigned i = 0; i < num_nodes; ++i)
+ {
+ // TODO PETSc negative indices?
+ auto const idx = num_components * global_indices[i] + comp;
+ nodal_vals_element[i] = (*_nodal_values)[idx];
+ }
- _residuals.set(element_index, std::sqrt(residual / num_int_pts));
+ double const residual = (interpolation_matrix * nodal_vals_element -
+ int_pt_vals_mat.row(comp).transpose())
+ .squaredNorm();
+
+ auto const eidx = num_components * element_index + comp;
+ _residuals->set(eidx, std::sqrt(residual / num_int_pts));
+ }
}
} // namespace NumLib
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index aa8929267d47411a80203c52aec28ee973a1013e..f1cbb79f266c7ef8ee355496cc4459ff62430be1 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
@@ -48,8 +48,11 @@ public:
explicit LocalLinearLeastSquaresExtrapolator(
NumLib::LocalToGlobalIndexMap const& dof_table);
- void extrapolate(
- ExtrapolatableElementCollection const& extrapolatables) override;
+ void extrapolate(const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table) override;
/*! \copydoc Extrapolator::calculateResiduals()
*
@@ -59,41 +62,44 @@ public:
* again.
*/
void calculateResiduals(
- ExtrapolatableElementCollection const& extrapolatables) override;
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table) override;
GlobalVector const& getNodalValues() const override
{
- return _nodal_values;
+ return *_nodal_values;
}
GlobalVector const& getElementResiduals() const override
{
- return _residuals;
- }
-
- ~LocalLinearLeastSquaresExtrapolator() override
- {
- NumLib::GlobalVectorProvider::provider.releaseVector(
- _nodal_values);
+ return *_residuals;
}
private:
//! Extrapolate one element.
void extrapolateElement(
- std::size_t const element_index,
- ExtrapolatableElementCollection const& extrapolatables,
- GlobalVector& counts);
+ std::size_t const element_index, const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables, const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table, GlobalVector& counts);
//! Compute the residuals for one element
void calculateResidualElement(
std::size_t const element_index,
- ExtrapolatableElementCollection const& extrapolatables);
+ const unsigned num_components,
+ ExtrapolatableElementCollection const& extrapolatables,
+ const double t,
+ GlobalVector const& current_solution,
+ LocalToGlobalIndexMap const& dof_table);
- GlobalVector& _nodal_values; //!< extrapolated nodal values
- GlobalVector _residuals; //!< extrapolation residuals
+ std::unique_ptr<GlobalVector> _nodal_values; //!< extrapolated nodal values
+ std::unique_ptr<GlobalVector> _residuals; //!< extrapolation residuals
//! DOF table used for writing to global vectors.
- NumLib::LocalToGlobalIndexMap const& _local_to_global;
+ NumLib::LocalToGlobalIndexMap const& _dof_table_single_component;
//! Avoids frequent reallocations.
std::vector<double> _integration_point_values_cache;
@@ -108,12 +114,12 @@ private:
Eigen::MatrixXd A_pinv;
};
- /*! Maps (\#nodes, \#int_pts) to (N_0, QR decomposition), where N_0 is the
- * shape matrix of the first integration point.
+ /*! Maps (#nodes, #int_pts) to (N_0, QR decomposition),
+ * where N_0 is the shape matrix of the first integration point.
*
- * \note It is assumed that the pair (\#nodes, \#int_pts) uniquely
- * identifies the set of all shape matrices N for a mesh element (i.e., only
- * N, not dN/dx etc.).
+ * \note It is assumed that the pair (#nodes, #int_pts) uniquely identifies
+ * the set of all shape matrices N for a mesh element (i.e., only N, not
+ * dN/dx etc.).
*
* \todo Add the element dimension as identifying criterion, or change to
* typeid.