diff --git a/NumLib/Extrapolation/ExtrapolatableElementCollection.h b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
index c07d719a79a71a8cfd4e3c7d1bb53f0312c9f582..e2d861360e913ed59e0bb9ad7b2e52e1b6ade995 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).
*
@@ -31,6 +37,11 @@ public:
/*! Returns integration point values of some property of a specific element.
*
* \param id ID of the element of which the property is queried.
+ * \param t The time used in the evaluation if time dependent quantities.
+ * \param current_solution The current solution of a ProcessLib::Process;
+ * more generally any nodal GlobalVector.
+ * \param dof_table The processes d.o.f. table used to get each element's
+ * local d.o.f. from \c current_solution.
* \param cache Can be used to compute a property on-the-fly.
*
* \remark
@@ -50,7 +61,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 +78,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 +94,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 +110,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 +124,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..3c483698556859d67b7366927f8a905a2e04f803 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,113 @@ 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);
- counts->setZero(); // TODO BLAS?
+ MathLib::MatrixVectorTraits<GlobalVector>::newInstance(*_nodal_values);
+ counts->setZero();
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();
+ auto const num_nodes = static_cast<unsigned>(N_0.cols());
+ auto const num_values =
+ static_cast<unsigned>(integration_point_values.size());
- assert(num_int_pts >= num_nodes &&
- "Least squares is not possible if there are more nodes than"
- "integration points.");
+ 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);
+
+ // 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{});
@@ -111,7 +155,8 @@ void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
interpolation_matrix.resize(num_int_pts, num_nodes);
interpolation_matrix.row(0) = N_0;
- for (unsigned int_pt = 1; int_pt < num_int_pts; ++int_pt) {
+ for (unsigned int_pt = 1; int_pt < num_int_pts; ++int_pt)
+ {
auto const& shp_mat =
extrapolatables.getShapeMatrix(element_index, int_pt);
assert(shp_mat.cols() == num_nodes);
@@ -122,7 +167,8 @@ void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
// JacobiSVD is extremely reliable, but fast only for small matrices.
// But we usually have small matrices and we don't compute very often.
- // Cf. http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
+ // Cf.
+ // http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
//
// Decomposes interpolation_matrix = U S V^T.
Eigen::JacobiSVD<Eigen::MatrixXd> svd(
@@ -137,60 +183,115 @@ void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
assert(rank == num_nodes);
// cf. http://eigen.tuxfamily.org/dox/JacobiSVD_8h_source.html
- cached_data.A_pinv.noalias() =
- V.leftCols(rank) *
- S.head(rank).asDiagonal().inverse() *
- U.leftCols(rank).transpose();
+ cached_data.A_pinv.noalias() = V.leftCols(rank) *
+ S.head(rank).asDiagonal().inverse() *
+ U.leftCols(rank).transpose();
}
- else if (cached_data.A.row(0) != N_0) {
+ else if (cached_data.A.row(0) != N_0)
+ {
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);
+ auto const& global_indices =
+ _dof_table_single_component(element_index, 0).rows;
- // 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;
+ if (num_components == 1)
+ {
+ 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? E.g., writing to ghost
+ // nodes? Furthermore: Is ghost nodes communication necessary for PETSc?
+ _nodal_values->add(global_indices, nodal_values);
+ counts.add(global_indices,
+ std::vector<double>(global_indices.size(), 1.0));
+ }
+ else
+ {
+ 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());
+
+ // _nodal_values is ordered location-wise
+ 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;
+ auto const num_values = static_cast<unsigned>(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]];
- }
+ const auto& global_indices =
+ _dof_table_single_component(element_index, 0).rows;
+ const auto num_nodes = static_cast<unsigned>(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);
- _residuals.set(element_index, std::sqrt(residual / 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];
+ }
+
+ double const residual = (interpolation_matrix * nodal_vals_element -
+ int_pt_vals_mat.row(comp).transpose())
+ .squaredNorm();
+
+ auto const eidx =
+ static_cast<GlobalIndexType>(num_components * element_index + comp);
+ // The residual is set to the root mean square value.
+ auto const root_mean_square = std::sqrt(residual / num_int_pts);
+ _residuals->set(eidx, root_mean_square);
+ }
}
} // namespace NumLib
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index aa8929267d47411a80203c52aec28ee973a1013e..11727b25e1d25ae60b1348ed85802c80cbe4b052 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
@@ -28,12 +28,10 @@ namespace NumLib
* the residuals are computed from the actual interpolation result that is being
* returned by this class.
*
- * Currently this class only supports interpolating single component variables.
- *
- * Furthermore, the number of integration points in each element must be greater
- * than or equal to the number of nodes of that element. This restriction is due to
- * the use of the least squares which requires an exact or overdetermined equation
- * system.
+ * The number of integration points in each element must be greater than or
+ * equal to the number of nodes of that element. This restriction is due
+ * to the use of the least squares which requires an exact or overdetermined
+ * equation system.
* \endparblock
*/
class LocalLinearLeastSquaresExtrapolator : public Extrapolator
@@ -48,8 +46,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 +60,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,8 +112,8 @@ 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
diff --git a/ProcessLib/CachedSecondaryVariable.cpp b/ProcessLib/CachedSecondaryVariable.cpp
index 3b29c6f9abb6907010db53fde690decb4ff81bb6..a5d442ca4111a98e4935f9629043268a578b8458 100644
--- a/ProcessLib/CachedSecondaryVariable.cpp
+++ b/ProcessLib/CachedSecondaryVariable.cpp
@@ -32,23 +32,27 @@ SecondaryVariableFunctions CachedSecondaryVariable::getExtrapolator()
{
// TODO copied from makeExtrapolator()
auto const eval_residuals = [this](
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ const double t,
+ GlobalVector const& x,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector> & /*result_cache*/
) -> GlobalVector const& {
- _extrapolator.calculateResiduals(*_extrapolatables);
+ _extrapolator.calculateResiduals(1, *_extrapolatables, t, x, dof_table);
return _extrapolator.getElementResiduals();
};
- return {BaseLib::easyBind(&CachedSecondaryVariable::evalField, this),
+ return {1, BaseLib::easyBind(&CachedSecondaryVariable::evalField, this),
eval_residuals};
}
GlobalVector const& CachedSecondaryVariable::evalField(
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ const double t,
+ GlobalVector const& x,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& /*result_cache*/
) const
{
+ (void)t, (void)x, (void)dof_table;
+ assert(t == _t && &x == _current_solution && &dof_table == _dof_table);
return evalFieldNoArgs();
}
@@ -60,7 +64,8 @@ GlobalVector const& CachedSecondaryVariable::evalFieldNoArgs() const
return _cached_nodal_values;
}
DBUG("Recomputing %s.", _internal_variable_name.c_str());
- _extrapolator.extrapolate(*_extrapolatables);
+ _extrapolator.extrapolate(
+ 1, *_extrapolatables, _t, *_current_solution, *_dof_table);
auto const& nodal_values = _extrapolator.getNodalValues();
MathLib::LinAlg::copy(nodal_values, _cached_nodal_values);
_needs_recomputation = false;
diff --git a/ProcessLib/CachedSecondaryVariable.h b/ProcessLib/CachedSecondaryVariable.h
index 0c6b5335f6e8c1353cf7cc7343582033baa71d47..28fe3a6e09bd81016cba33891ab2b2c36cb8b99b 100644
--- a/ProcessLib/CachedSecondaryVariable.h
+++ b/ProcessLib/CachedSecondaryVariable.h
@@ -60,8 +60,19 @@ public:
//! Returns extrapolation functions that compute the secondary variable.
SecondaryVariableFunctions getExtrapolator();
- //! Set that recomputation is necessary.
- void expire() { _needs_recomputation = true; }
+ void setTime(const double t)
+ {
+ _t = t;
+ _needs_recomputation = true;
+ }
+
+ void updateCurrentSolution(GlobalVector const& x,
+ NumLib::LocalToGlobalIndexMap const& dof_table)
+ {
+ _current_solution = &x;
+ _dof_table = &dof_table;
+ _needs_recomputation = true;
+ }
private:
//! Provides the value at the current index of the _context.
@@ -69,8 +80,9 @@ private:
//! Computes the secondary Variable.
GlobalVector const& evalField(
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ const double t,
+ GlobalVector const& x,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& /*result_cache*/
) const;
@@ -85,6 +97,10 @@ private:
std::unique_ptr<NumLib::ExtrapolatableElementCollection> _extrapolatables;
SecondaryVariableContext const& _context;
std::string const _internal_variable_name;
+
+ double _t = 0.0;
+ GlobalVector const* _current_solution = nullptr;
+ NumLib::LocalToGlobalIndexMap const* _dof_table = nullptr;
};
} // namespace ProcessLib
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 852a8efbaaf547639309e5b1d53345c00fd7a730..09ea5a4bc103de7cc874b77f23ea964c1d119d81 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -50,12 +50,21 @@ class ComponentTransportLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -320,6 +329,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 0);
@@ -327,6 +339,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 1);
@@ -334,6 +349,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 2);
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index 3b83bc3d4e405d5d5d5ab9704c705c40d336d6e1..2ebabe39ebe8fa8b95d86bf10dd2f08864ea0a2d 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -45,24 +45,24 @@ void ComponentTransportProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "darcy_velocity_x",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ComponentTransportLocalAssemblerInterface::
getIntPtDarcyVelocityX));
if (mesh.getDimension() > 1)
{
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "darcy_velocity_y",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ComponentTransportLocalAssemblerInterface::
getIntPtDarcyVelocityY));
}
if (mesh.getDimension() > 2)
{
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "darcy_velocity_z",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ComponentTransportLocalAssemblerInterface::
getIntPtDarcyVelocityZ));
}
diff --git a/ProcessLib/GlobalVectorFromNamedFunction.cpp b/ProcessLib/GlobalVectorFromNamedFunction.cpp
index 82020a384e3361a22acf8eed09bd3cee2ee338d2..9a456d918e547dffdf0632d2bca74337d6c3ef16 100644
--- a/ProcessLib/GlobalVectorFromNamedFunction.cpp
+++ b/ProcessLib/GlobalVectorFromNamedFunction.cpp
@@ -28,6 +28,7 @@ GlobalVectorFromNamedFunction::GlobalVectorFromNamedFunction(
}
GlobalVector const& GlobalVectorFromNamedFunction::call(
+ const double /*t*/,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result)
diff --git a/ProcessLib/GlobalVectorFromNamedFunction.h b/ProcessLib/GlobalVectorFromNamedFunction.h
index 83ab3aa8e478d3d3e457c4e5c2298c8ffa899628..5b32bb993db64fe97dbba421eb10bcefc14f41b8 100644
--- a/ProcessLib/GlobalVectorFromNamedFunction.h
+++ b/ProcessLib/GlobalVectorFromNamedFunction.h
@@ -42,7 +42,7 @@ public:
//! The signature of this method matches
//! SecondaryVariableFunctions::Function, i.e., this method can be used to
//! compute a secondary variable.
- GlobalVector const& call(GlobalVector const& x,
+ GlobalVector const& call(const double t, GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result);
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index d9734be0888bb5b56d5a0f9a8be1b3dd2ff660d2..55f965a7d7b8283c2e7597ae66a15a258b71adcc 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -13,6 +13,7 @@
#include "GroundwaterFlowProcessData.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -32,13 +33,10 @@ class GroundwaterFlowLocalAssemblerInterface
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(
+ virtual std::vector<double> const& getIntPtDarcyVelocity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -69,10 +67,7 @@ public:
_integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod, GlobalDim>(
- element, is_axially_symmetric, _integration_method)),
- _darcy_velocities(
- GlobalDim,
- std::vector<double>(_integration_method.getNumberOfPoints()))
+ element, is_axially_symmetric, _integration_method))
{
}
@@ -107,39 +102,6 @@ public:
}
}
- void computeSecondaryVariableConcrete(
- const double t,
- std::vector<double> const& local_x) override
- {
- auto const local_matrix_size = local_x.size();
- // This assertion is valid only if all nodal d.o.f. use the same shape
- // matrices.
- assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- const auto local_x_vec =
- MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
-
- SpatialPosition pos;
- pos.setElementID(_element.getID());
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- pos.setIntegrationPoint(ip);
- auto const& sm = _shape_matrices[ip];
- auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
- // Darcy velocity only computed for output.
- GlobalDimVectorType const darcy_velocity = -k * sm.dNdx * local_x_vec;
-
- for (unsigned d = 0; d < GlobalDim; ++d)
- {
- _darcy_velocities[d][ip] = darcy_velocity[d];
- }
- }
- }
-
/// Computes the flux in the point \c p_local_coords that is given in local
/// coordinates using the values from \c local_x.
// TODO add time dependency
@@ -187,25 +149,40 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
- std::vector<double> const& getIntPtDarcyVelocityX(
- std::vector<double>& /*cache*/) const override
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const override
{
- assert(!_darcy_velocities.empty());
- return _darcy_velocities[0];
- }
+ auto const n_integration_points =
+ _integration_method.getNumberOfPoints();
- std::vector<double> const& getIntPtDarcyVelocityY(
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 1);
- return _darcy_velocities[1];
- }
+ auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto const local_x = current_solution.get(indices);
+ auto const local_x_vec =
+ MathLib::toVector<Eigen::Matrix<double, ShapeFunction::NPOINTS, 1>>(
+ local_x, ShapeFunction::NPOINTS);
- std::vector<double> const& getIntPtDarcyVelocityZ(
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 2);
- return _darcy_velocities[2];
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, GlobalDim, n_integration_points);
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+
+ for (unsigned i = 0; i < n_integration_points; ++i)
+ {
+ pos.setIntegrationPoint(i);
+ auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
+ // dimensions: (d x 1) = (d x n) * (n x 1)
+ cache_mat.col(i).noalias() =
+ -k * _shape_matrices[i].dNdx * local_x_vec;
+ }
+
+ return cache;
}
private:
@@ -215,8 +192,6 @@ private:
IntegrationMethod const _integration_method;
std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
_shape_matrices;
-
- std::vector<std::vector<double>> _darcy_velocities;
};
} // namespace GroundwaterFlow
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index d50fb4e08c6acaeaca331d956be82fff92f83660..010552376497d308aa31f8877c2c94fc65f25a87 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -51,25 +51,10 @@ void GroundwaterFlowProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
+ "darcy_velocity",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
- &GroundwaterFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
- if (mesh.getDimension() > 1) {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &GroundwaterFlowLocalAssemblerInterface::
- getIntPtDarcyVelocityY));
- }
- if (mesh.getDimension() > 2) {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &GroundwaterFlowLocalAssemblerInterface::
- getIntPtDarcyVelocityZ));
- }
+ mesh.getDimension(), getExtrapolator(), _local_assemblers,
+ &GroundwaterFlowLocalAssemblerInterface::getIntPtDarcyVelocity));
}
void GroundwaterFlowProcess::assembleConcreteProcess(const double t,
@@ -103,16 +88,5 @@ void GroundwaterFlowProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac, coupling_term);
}
-
-void GroundwaterFlowProcess::computeSecondaryVariableConcrete(
- const double t, GlobalVector const& x,
- StaggeredCouplingTerm const& coupled_term)
-{
- DBUG("Compute the velocity for GroundwaterFlowProcess.");
- GlobalExecutor::executeMemberOnDereferenced(
- &GroundwaterFlowLocalAssemblerInterface::computeSecondaryVariable,
- _local_assemblers, *_local_to_global_index_map, t, x, coupled_term);
-}
-
} // namespace GroundwaterFlow
} // namespace ProcessLib
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index eadd1739792ce1aa9abfb4827f5ce130d743e213..43a783e95ebd3500e9d693bb829434fb057ae0cd 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -87,11 +87,6 @@ public:
}
}
- void computeSecondaryVariableConcrete(double const t,
- GlobalVector const& x,
- StaggeredCouplingTerm const&
- coupled_term) override;
-
private:
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 6eca1ffe560a61310356f7605e35abb3113addaa..bf928cb0a57a4081285fd2829a3f543359fc0235 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -15,6 +15,7 @@
#include "HTProcessData.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -50,13 +51,10 @@ class HTLocalAssemblerInterface
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(
+ virtual std::vector<double> const& getIntPtDarcyVelocity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -90,10 +88,7 @@ public:
HTProcessData const& process_data)
: _element(element),
_process_data(process_data),
- _integration_method(integration_order),
- _darcy_velocities(
- GlobalDim,
- std::vector<double>(_integration_method.getNumberOfPoints()))
+ _integration_method(integration_order)
{
// This assertion is valid only if all nodal d.o.f. use the same shape
// matrices.
@@ -263,20 +258,38 @@ public:
}
}
- void computeSecondaryVariableConcrete(
- double const t, std::vector<double> const& local_x) override
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _ip_data[integration_point].N;
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double t,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const override
{
+ auto const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto const local_x = current_solution.get(indices);
+
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, GlobalDim, n_integration_points);
+
SpatialPosition pos;
pos.setElementID(_element.getID());
- auto const K =
- _process_data.porous_media_properties.getIntrinsicPermeability(t,
- pos);
MaterialLib::Fluid::FluidProperty::ArrayType vars;
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
@@ -294,55 +307,27 @@ public:
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
+ auto const K =
+ _process_data.porous_media_properties.getIntrinsicPermeability(
+ t, pos);
+
auto const mu = _process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
GlobalDimMatrixType const K_over_mu = K / mu;
- GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values;
+ cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
+
if (_process_data.has_gravity)
{
auto const rho_w = _process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Density, vars);
auto const b = _process_data.specific_body_force;
// here it is assumed that the vector b is directed 'downwards'
- velocity += K_over_mu * rho_w * b;
- }
-
- for (unsigned d = 0; d < GlobalDim; ++d)
- {
- _darcy_velocities[d][ip] = velocity[d];
+ cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
}
}
- }
-
- Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
- const unsigned integration_point) const override
- {
- auto const& N = _ip_data[integration_point].N;
- // assumes N is stored contiguously in memory
- return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
- }
-
- std::vector<double> const& getIntPtDarcyVelocityX(
- std::vector<double>& /*cache*/) const override
- {
- assert(!_darcy_velocities.empty());
- return _darcy_velocities[0];
- }
-
- std::vector<double> const& getIntPtDarcyVelocityY(
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 1);
- return _darcy_velocities[1];
- }
-
- std::vector<double> const& getIntPtDarcyVelocityZ(
- std::vector<double>& /*cache*/) const override
- {
- assert(_darcy_velocities.size() > 2);
- return _darcy_velocities[2];
+ return cache;
}
private:
@@ -355,7 +340,6 @@ private:
Eigen::aligned_allocator<
IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
_ip_data;
- std::vector<std::vector<double>> _darcy_velocities;
};
} // namespace HT
diff --git a/ProcessLib/HT/HTProcess.cpp b/ProcessLib/HT/HTProcess.cpp
index 7a70c5ad7671133cc2b60e663432a0bba279e2f8..8da71f612fdfd3e377fa69c2e66face6d23af661 100644
--- a/ProcessLib/HT/HTProcess.cpp
+++ b/ProcessLib/HT/HTProcess.cpp
@@ -45,26 +45,10 @@ void HTProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &HTLocalAssemblerInterface::getIntPtDarcyVelocityX));
-
- if (mesh.getDimension() > 1)
- {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
- makeExtrapolator(
- getExtrapolator(), _local_assemblers,
- &HTLocalAssemblerInterface::getIntPtDarcyVelocityY));
- }
- if (mesh.getDimension() > 2)
- {
- _secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
- makeExtrapolator(
- getExtrapolator(), _local_assemblers,
- &HTLocalAssemblerInterface::getIntPtDarcyVelocityZ));
- }
+ "darcy_velocity",
+ makeExtrapolator(mesh.getDimension(), getExtrapolator(),
+ _local_assemblers,
+ &HTLocalAssemblerInterface::getIntPtDarcyVelocity));
}
void HTProcess::assembleConcreteProcess(const double t,
@@ -98,16 +82,6 @@ void HTProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac, coupling_term);
}
-void HTProcess::computeSecondaryVariableConcrete(
- double const t, GlobalVector const& x,
- StaggeredCouplingTerm const& coupling_term)
-{
- DBUG("Compute the Darcy velocity for HTProcess.");
- GlobalExecutor::executeMemberOnDereferenced(
- &HTLocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
- *_local_to_global_index_map, t, x, coupling_term);
-}
-
} // namespace HT
} // namespace ProcessLib
diff --git a/ProcessLib/HT/HTProcess.h b/ProcessLib/HT/HTProcess.h
index 6dcdaaeaf2669de46918240013aba189f7687cdb..a5048e295cd1daaa2f4507d3d955a06bb0a8d8d6 100644
--- a/ProcessLib/HT/HTProcess.h
+++ b/ProcessLib/HT/HTProcess.h
@@ -56,9 +56,6 @@ public:
bool isLinear() const override { return false; }
//! @}
- void computeSecondaryVariableConcrete(
- double const t, GlobalVector const& x,
- StaggeredCouplingTerm const& coupling_term) override;
private:
void initializeConcreteProcess(
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index 1dc74f69749763d78eee31d073304c1e49ee5b1e..4c2070c6d65d091a106f041d7fce6bd9b64589b0 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -37,12 +37,21 @@ class HeatConductionLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtHeatFluxX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtHeatFluxY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtHeatFluxZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -169,6 +178,9 @@ public:
}
std::vector<double> const& getIntPtHeatFluxX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_heat_fluxes.empty());
@@ -176,6 +188,9 @@ public:
}
std::vector<double> const& getIntPtHeatFluxY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_heat_fluxes.size() > 1);
@@ -183,6 +198,9 @@ public:
}
std::vector<double> const& getIntPtHeatFluxZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_heat_fluxes.size() > 2);
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index d341c30a4f9ab2476aa4209820f4b591b3e90821..b2e8e3a47422167747d6a232aadab614cd5c61d1 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -61,25 +61,25 @@ void HeatConductionProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "heat_flux_x", 1,
+ "heat_flux_x",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HeatConductionLocalAssemblerInterface::getIntPtHeatFluxX));
if (mesh.getDimension() > 1)
{
_secondary_variables.addSecondaryVariable(
- "heat_flux_y", 1,
+ "heat_flux_y",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HeatConductionLocalAssemblerInterface::getIntPtHeatFluxY));
}
if (mesh.getDimension() > 2)
{
_secondary_variables.addSecondaryVariable(
- "heat_flux_z", 1,
+ "heat_flux_z",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HeatConductionLocalAssemblerInterface::getIntPtHeatFluxZ));
}
}
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index fb84183887b128d3dab10a69c2024a378ceb293f..ad1c78ca441630a0ee266bb661919d75841e3e85 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -102,48 +102,93 @@ struct HydroMechanicsLocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
virtual std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
};
@@ -462,30 +507,45 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 1);
}
std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 3);
}
std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -493,6 +553,9 @@ public:
}
std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -500,30 +563,45 @@ public:
}
std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 0);
}
std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 1);
}
std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 2);
}
std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 3);
}
std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -531,6 +609,9 @@ public:
}
std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -538,6 +619,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_darcy_velocities.empty());
@@ -545,6 +629,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 1);
@@ -552,6 +639,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 2);
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index e4ae10e2839e15d19b03003a95038e39ba11c627..5211eb2345372fbacaba405b4c67931bb136edda 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -114,83 +114,83 @@ private:
NumLib::ComponentOrder::BY_LOCATION);
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xx", 1,
+ "sigma_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtSigmaXX));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yy", 1,
+ "sigma_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtSigmaYY));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_zz", 1,
+ "sigma_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtSigmaZZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xy", 1,
+ "sigma_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtSigmaXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xz", 1,
+ "sigma_xz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtSigmaXZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yz", 1,
+ "sigma_yz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtSigmaYZ));
}
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xx", 1,
+ "epsilon_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtEpsilonXX));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yy", 1,
+ "epsilon_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtEpsilonYY));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_zz", 1,
+ "epsilon_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtEpsilonZZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xy", 1,
+ "epsilon_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::getIntPtEpsilonXY));
Base::_secondary_variables.addSecondaryVariable(
- "velocity_x", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "velocity_x",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::
getIntPtDarcyVelocityX));
Base::_secondary_variables.addSecondaryVariable(
- "velocity_y", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "velocity_y",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::
getIntPtDarcyVelocityY));
Base::_secondary_variables.addSecondaryVariable(
- "velocity_z", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "velocity_z",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&HydroMechanicsLocalAssemblerInterface::
getIntPtDarcyVelocityZ));
}
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h
index 928b52e286185eb4fce305c7cc9dcf44f2a26c96..6a4a398b18bfd4f72bac6e3b3581c82b1cac55e2 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h
@@ -92,7 +92,6 @@ public:
void postTimestepConcrete(std::vector<double> const& /*local_x*/) override;
-
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -103,30 +102,45 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 1);
}
std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 3);
}
std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -134,6 +148,9 @@ public:
}
std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -141,30 +158,45 @@ public:
}
std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 0);
}
std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 1);
}
std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 2);
}
std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 3);
}
std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -172,6 +204,9 @@ public:
}
std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerInterface.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerInterface.h
index 3cf4f7410d18af855e21e62ea94c08505a25cd79..907b70fa56cfb215c81327866c8b6be7a16d4faf 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerInterface.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerInterface.h
@@ -83,39 +83,75 @@ public:
}
virtual std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
private:
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
index e9cbc96638813d20d1f04ca13f03b59ae63004f8..b8eb7ff4cb96acb08ba67cb87bb3b0f279c63936 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
@@ -102,30 +102,45 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 1);
}
std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 3);
}
std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -133,6 +148,9 @@ public:
}
std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -140,30 +158,45 @@ public:
}
std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 0);
}
std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 1);
}
std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 2);
}
std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 3);
}
std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -171,6 +204,9 @@ public:
}
std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
index 50c8b4b7ff2e5a2e27952433077658bda53fd8f9..a91d872156526a94ac1bce973df5338d26aba54d 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
@@ -104,30 +104,45 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 1);
}
std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 3);
}
std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -135,6 +150,9 @@ public:
}
std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -142,30 +160,45 @@ public:
}
std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 0);
}
std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 1);
}
std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 2);
}
std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 3);
}
std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -173,6 +206,9 @@ public:
}
std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index 89854b68d2f4652df2a3ae908d3b13949b1263d4..1035c8e2214c37a8bd46be6d9918a82dec07701b 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -186,80 +186,80 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
NumLib::ComponentOrder::BY_LOCATION);
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xx", 1,
+ "sigma_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaXX));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yy", 1,
+ "sigma_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaYY));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_zz", 1,
+ "sigma_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaZZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xy", 1,
+ "sigma_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xz", 1,
+ "sigma_xz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaXZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yz", 1,
+ "sigma_yz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaYZ));
}
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xx", 1,
+ "epsilon_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXX));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yy", 1,
+ "epsilon_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonYY));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_zz", 1,
+ "epsilon_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonZZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xy", 1,
+ "epsilon_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xz", 1,
+ "epsilon_xz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yz", 1,
+ "epsilon_yz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonYZ));
}
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index c6adbf2bf321e66780af45afb02accca38680ec8..b2c89ddc270d28aa48e2b0b5195f5f7ba0b5160f 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -40,12 +40,21 @@ class LiquidFlowLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -114,6 +123,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_darcy_velocities.empty());
@@ -121,6 +133,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 1);
@@ -128,6 +143,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 2);
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index 3eaaa0b739b2a0db75106c0349013255427699a8..6b61dc812a5721c3da77be9820e4b38907bca094 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -69,25 +69,25 @@ void LiquidFlowProcess::initializeConcreteProcess(
*_material_properties);
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
+ "darcy_velocity_x",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
if (mesh.getDimension() > 1)
{
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
+ "darcy_velocity_y",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
}
if (mesh.getDimension() > 2)
{
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
+ "darcy_velocity_z",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
}
}
diff --git a/ProcessLib/Output.cpp b/ProcessLib/Output.cpp
index 2e2cc19a789521a636b524a94bc908a635efd178..969abda66d2523e918131576a2bbcbb13e0e4190 100644
--- a/ProcessLib/Output.cpp
+++ b/ProcessLib/Output.cpp
@@ -121,7 +121,7 @@ void Output::doOutputAlways(Process const& process,
+ ".vtu";
std::string const output_file_path = BaseLib::joinPaths(_output_directory, output_file_name);
DBUG("output to %s", output_file_path.c_str());
- doProcessOutput(output_file_path, _output_file_compression, x,
+ doProcessOutput(output_file_path, _output_file_compression, t, x,
process.getMesh(), process.getDOFTable(),
process.getProcessVariables(),
process.getSecondaryVariables(), process_output);
@@ -185,7 +185,7 @@ void Output::doOutputNonlinearIteration(Process const& process,
+ ".vtu";
std::string const output_file_path = BaseLib::joinPaths(_output_directory, output_file_name);
DBUG("output iteration results to %s", output_file_path.c_str());
- doProcessOutput(output_file_path, _output_file_compression, x,
+ doProcessOutput(output_file_path, _output_file_compression, t, x,
process.getMesh(), process.getDOFTable(),
process.getProcessVariables(),
process.getSecondaryVariables(), process_output);
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 25e265394bbc0180413bf33b88b40bec6bbc273e..81a7084b1fabde10cd2217c09df70b67dab03d4c 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -217,16 +217,14 @@ void Process::finishNamedFunctionsInitialization()
for (auto const& named_function :
_named_function_caller.getNamedFunctions()) {
auto const& name = named_function.getName();
- // secondary variables generated from named functions have the prefix
- // "fct_".
_secondary_variables.addSecondaryVariable(
- "fct_" + name, 1,
- {BaseLib::easyBind(
- &GlobalVectorFromNamedFunction::call,
- GlobalVectorFromNamedFunction(
- _named_function_caller.getSpecificFunctionCaller(name), _mesh,
- getSingleComponentDOFTable(),
- _secondary_variable_context)),
+ name,
+ {1, BaseLib::easyBind(
+ &GlobalVectorFromNamedFunction::call,
+ GlobalVectorFromNamedFunction(
+ _named_function_caller.getSpecificFunctionCaller(name),
+ _mesh, getSingleComponentDOFTable(),
+ _secondary_variable_context)),
nullptr});
}
}
@@ -240,6 +238,11 @@ void Process::computeSparsityPattern()
void Process::preTimestep(GlobalVector const& x, const double t,
const double delta_t)
{
+ for (auto& cached_var : _cached_secondary_variables)
+ {
+ cached_var->setTime(t);
+ }
+
MathLib::LinAlg::setLocalAccessibleVector(x);
preTimestepConcreteProcess(x, t, delta_t);
_boundary_conditions.preTimestep(t);
@@ -264,7 +267,7 @@ void Process::preIteration(const unsigned iter, const GlobalVector &x)
{
// In every new iteration cached values of secondary variables are expired.
for (auto& cached_var : _cached_secondary_variables) {
- cached_var->expire();
+ cached_var->updateCurrentSolution(x, *_local_to_global_index_map);
}
MathLib::LinAlg::setLocalAccessibleVector(x);
diff --git a/ProcessLib/ProcessOutput.cpp b/ProcessLib/ProcessOutput.cpp
index 7f90dcacb1128f79b6344cf2da1e33786294ed8f..5daa20d6aa3d340461a949620b57ef69e12173b7 100644
--- a/ProcessLib/ProcessOutput.cpp
+++ b/ProcessLib/ProcessOutput.cpp
@@ -42,6 +42,7 @@ ProcessOutput::ProcessOutput(BaseLib::ConfigTree const& output_config)
void doProcessOutput(std::string const& file_name,
bool const compress_output,
+ const double t,
GlobalVector const& x,
MeshLib::Mesh& mesh,
NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -117,27 +118,43 @@ void doProcessOutput(std::string const& file_name,
}
#ifndef USE_PETSC
- // the following section is for the output of secondary variables
-
auto add_secondary_var = [&](SecondaryVariable const& var,
- std::string const& output_name)
- {
- assert(var.n_components == 1); // TODO implement other cases
-
+ std::string const& output_name) {
{
DBUG(" secondary variable %s", output_name.c_str());
- auto result = MeshLib::getOrCreateMeshProperty<double>(
- mesh, output_name, MeshLib::MeshItemType::Node, 1);
+ auto& nodal_values_mesh = *MeshLib::getOrCreateMeshProperty<double>(
+ mesh, output_name, MeshLib::MeshItemType::Node,
+ var.fcts.num_components);
+ if (nodal_values_mesh.size() !=
+ mesh.getNumberOfNodes() * var.fcts.num_components)
+ {
+ OGS_FATAL(
+ "Nodal property `%s' does not have the right number of "
+ "components. Expected: %d, actual: %d",
+ output_name.c_str(),
+ mesh.getNumberOfNodes() * var.fcts.num_components,
+ nodal_values_mesh.size());
+ }
std::unique_ptr<GlobalVector> result_cache;
auto const& nodal_values =
- var.fcts.eval_field(x, dof_table, result_cache);
+ var.fcts.eval_field(t, x, dof_table, result_cache);
+ if (nodal_values_mesh.size() !=
+ static_cast<std::size_t>(nodal_values.size()))
+ {
+ OGS_FATAL(
+ "Secondary variable `%s' did not evaluate to the right "
+ "number of components. Expected: %d, actual: %d.",
+ var.name.c_str(), nodal_values_mesh.size(),
+ nodal_values.size());
+ }
// Copy result
- for (GlobalIndexType i = 0; i < nodal_values.size(); ++i) {
+ for (GlobalIndexType i = 0; i < nodal_values.size(); ++i)
+ {
assert(!std::isnan(nodal_values[i]));
- (*result)[i] = nodal_values[i];
+ nodal_values_mesh[i] = nodal_values[i];
}
}
@@ -146,18 +163,38 @@ void doProcessOutput(std::string const& file_name,
DBUG(" secondary variable %s residual", output_name.c_str());
auto const& property_name_res = output_name + "_residual";
- auto result = MeshLib::getOrCreateMeshProperty<double>(
- mesh, property_name_res, MeshLib::MeshItemType::Cell, 1);
+ auto& residuals_mesh = *MeshLib::getOrCreateMeshProperty<double>(
+ mesh, property_name_res, MeshLib::MeshItemType::Cell,
+ var.fcts.num_components);
+ if (residuals_mesh.size() !=
+ mesh.getNumberOfElements() * var.fcts.num_components)
+ {
+ OGS_FATAL(
+ "Cell property `%s' does not have the right number of "
+ "components. Expected: %d, actual: %d",
+ property_name_res.c_str(),
+ mesh.getNumberOfElements() * var.fcts.num_components,
+ residuals_mesh.size());
+ }
std::unique_ptr<GlobalVector> result_cache;
auto const& residuals =
- var.fcts.eval_residuals(x, dof_table, result_cache);
+ var.fcts.eval_residuals(t, x, dof_table, result_cache);
+ if (residuals_mesh.size() !=
+ static_cast<std::size_t>(residuals.size()))
+ {
+ OGS_FATAL(
+ "Thee residual of secondary variable `%s' did not evaluate "
+ "to the right number of components. Expected: %d, actual: "
+ "%d.",
+ var.name.c_str(), residuals_mesh.size(), residuals.size());
+ }
// Copy result
for (GlobalIndexType i = 0; i < residuals.size(); ++i)
{
assert(!std::isnan(residuals[i]));
- (*result)[i] = residuals[i];
+ residuals_mesh[i] = residuals[i];
}
}
};
@@ -172,8 +209,6 @@ void doProcessOutput(std::string const& file_name,
add_secondary_var(secondary_variables.get(external_variable_name),
external_variable_name);
}
-
- // secondary variables output end
#else
(void) secondary_variables;
#endif // USE_PETSC
diff --git a/ProcessLib/ProcessOutput.h b/ProcessLib/ProcessOutput.h
index a1f5ff3c39af426acc294d3152f14abb6e46b601..e1e983c7e1d6a7d672f55dc2f20513a099a27a22 100644
--- a/ProcessLib/ProcessOutput.h
+++ b/ProcessLib/ProcessOutput.h
@@ -32,6 +32,7 @@ struct ProcessOutput final
//! Writes output to the given \c file_name using the VTU file format.
void doProcessOutput(std::string const& file_name,
bool const compress_output,
+ const double t,
GlobalVector const& x,
MeshLib::Mesh& mesh,
NumLib::LocalToGlobalIndexMap const& dof_table,
diff --git a/ProcessLib/RichardsFlow/RichardsFlowFEM.h b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
index 319bbacd977d968385ed72d5bb0d855934450262..08c61797ef38114dcd1d2a08e3d94286a6752baa 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowFEM.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
@@ -56,15 +56,27 @@ class RichardsFlowLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -294,6 +306,9 @@ public:
}
std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_saturation.empty());
@@ -301,6 +316,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_darcy_velocities.empty());
@@ -308,6 +326,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 1);
@@ -315,6 +336,9 @@ public:
}
std::vector<double> const& getIntPtDarcyVelocityZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 2);
diff --git a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
index c87d3dd6cdf843931c0a104d42f4ce22b5e3b5e8..073e95c291a64d0e7195a994287953c206d2719e 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
+++ b/ProcessLib/RichardsFlow/RichardsFlowProcess.cpp
@@ -49,31 +49,31 @@ void RichardsFlowProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "saturation", 1,
+ "saturation",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&RichardsFlowLocalAssemblerInterface::getIntPtSaturation));
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_x", 1,
+ "darcy_velocity_x",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
if (mesh.getDimension() > 1)
{
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_y", 1,
+ "darcy_velocity_y",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
}
if (mesh.getDimension() > 2)
{
_secondary_variables.addSecondaryVariable(
- "darcy_velocity_z", 1,
+ "darcy_velocity_z",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&RichardsFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
}
}
diff --git a/ProcessLib/SecondaryVariable.cpp b/ProcessLib/SecondaryVariable.cpp
index d25ef2032208c238373693a96fd0672767202c4e..ccbca39e5bced29610c9ddc85941530cb2c7c478 100644
--- a/ProcessLib/SecondaryVariable.cpp
+++ b/ProcessLib/SecondaryVariable.cpp
@@ -25,16 +25,15 @@ void SecondaryVariableCollection::addNameMapping(std::string const& internal_nam
}
void SecondaryVariableCollection::addSecondaryVariable(
- std::string const& internal_name,
- const unsigned num_components,
- SecondaryVariableFunctions&& fcts)
+ std::string const& internal_name, SecondaryVariableFunctions&& fcts)
{
if (!_configured_secondary_variables
.emplace(std::make_pair(
internal_name,
SecondaryVariable{internal_name /* TODO change */,
- num_components, std::move(fcts)}))
- .second) {
+ std::move(fcts)}))
+ .second)
+ {
OGS_FATAL(
"The secondary variable with internal name `%s' has already been "
"set up.",
diff --git a/ProcessLib/SecondaryVariable.h b/ProcessLib/SecondaryVariable.h
index de54dec8710a2c0378523a65a3ec576553e3ce02..11adea0e5b1aa8918f3aefa0312a5ea3a1e78396 100644
--- a/ProcessLib/SecondaryVariable.h
+++ b/ProcessLib/SecondaryVariable.h
@@ -37,58 +37,73 @@ struct SecondaryVariableFunctions final
* is stored somewhere else.
*/
using Function = std::function<GlobalVector const&(
+ const double t,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache)>;
SecondaryVariableFunctions() = default;
- template<typename F1, typename F2>
- SecondaryVariableFunctions(F1&& eval_field_, F2&& eval_residuals_)
- : eval_field(std::forward<F1>(eval_field_))
- , eval_residuals(std::forward<F2>(eval_residuals_))
+ template <typename F1, typename F2>
+ SecondaryVariableFunctions(const unsigned num_components_, F1&& eval_field_,
+ F2&& eval_residuals_)
+ : num_components(num_components_),
+ eval_field(std::forward<F1>(eval_field_)),
+ eval_residuals(std::forward<F2>(eval_residuals_))
{
// Used to detect nasty implicit conversions.
- static_assert(std::is_same<GlobalVector const&,
- typename std::result_of<F1(
- GlobalVector const&, NumLib::LocalToGlobalIndexMap const&,
- std::unique_ptr<GlobalVector>&
- )>::type>::value,
+ static_assert(
+ std::is_same<GlobalVector const&,
+ typename std::result_of<F1(
+ double const, GlobalVector const&,
+ NumLib::LocalToGlobalIndexMap const&,
+ std::unique_ptr<GlobalVector>&)>::type>::value,
"The function eval_field_ does not return a const reference"
" to a GlobalVector");
- static_assert(std::is_same<GlobalVector const&,
- typename std::result_of<F2(
- GlobalVector const&, NumLib::LocalToGlobalIndexMap const&,
- std::unique_ptr<GlobalVector>&
- )>::type>::value,
+ static_assert(
+ std::is_same<GlobalVector const&,
+ typename std::result_of<F2(
+ double const, GlobalVector const&,
+ NumLib::LocalToGlobalIndexMap const&,
+ std::unique_ptr<GlobalVector>&)>::type>::value,
"The function eval_residuals_ does not return a const reference"
" to a GlobalVector");
}
- template<typename F1>
- SecondaryVariableFunctions(F1&& eval_field_, std::nullptr_t)
- : eval_field(std::forward<F1>(eval_field_))
+ template <typename F1>
+ SecondaryVariableFunctions(const unsigned num_components_, F1&& eval_field_,
+ std::nullptr_t)
+ : num_components(num_components_),
+ eval_field(std::forward<F1>(eval_field_))
{
// Used to detect nasty implicit conversions.
- static_assert(std::is_same<GlobalVector const&,
- typename std::result_of<F1(
- GlobalVector const&, NumLib::LocalToGlobalIndexMap const&,
- std::unique_ptr<GlobalVector>&
- )>::type>::value,
+ static_assert(
+ std::is_same<GlobalVector const&,
+ typename std::result_of<F1(
+ double const, GlobalVector const&,
+ NumLib::LocalToGlobalIndexMap const&,
+ std::unique_ptr<GlobalVector>&)>::type>::value,
"The function eval_field_ does not return a const reference"
" to a GlobalVector");
}
- Function eval_field;
- Function eval_residuals;
+ const unsigned num_components; //!< Number of components of the variable.
+
+ //! Computes the value of the field at every node of the underlying mesh.
+ Function const eval_field;
+
+ //! If the secondary variable is extrapolated from integration points to
+ //! mesh nodes, this function computes the extrapolation residual. For
+ //! further information check the specific NumLib::Extrapolator
+ //! documentation.
+ Function const eval_residuals;
};
//! Stores information about a specific secondary variable
struct SecondaryVariable final
{
std::string const name; //!< Name of the variable; used, e.g., for output.
- const unsigned n_components; //!< Number of components of the variable.
//! Functions used for computing the secondary variable.
SecondaryVariableFunctions fcts;
@@ -106,7 +121,6 @@ public:
*
* \param internal_name the tag in the project file associated with this
* secondary variable.
- * \param num_components the variable's number of components.
* \param fcts functions that compute the variable.
*
* \note
@@ -115,7 +129,6 @@ public:
* All other variables are silently ignored.
*/
void addSecondaryVariable(std::string const& internal_name,
- const unsigned num_components,
SecondaryVariableFunctions&& fcts);
//! Returns the secondary variable with the given external name.
@@ -138,6 +151,7 @@ private:
/*! Creates an object that computes a secondary variable via extrapolation of
* integration point values.
*
+ * \param num_components The number of components of the secondary variable.
* \param extrapolator The extrapolator used for extrapolation.
* \param local_assemblers The collection of local assemblers whose integration
* point values will be extrapolated.
@@ -147,36 +161,42 @@ private:
*/
template <typename LocalAssemblerCollection>
SecondaryVariableFunctions makeExtrapolator(
+ const unsigned num_components,
NumLib::Extrapolator& extrapolator,
LocalAssemblerCollection const& local_assemblers,
typename NumLib::ExtrapolatableLocalAssemblerCollection<
LocalAssemblerCollection>::IntegrationPointValuesMethod
integration_point_values_method)
{
- auto const eval_field = [&extrapolator, &local_assemblers,
+ auto const eval_field = [num_components, &extrapolator, &local_assemblers,
integration_point_values_method](
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ const double t,
+ 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);
+ extrapolator.extrapolate(num_components, extrapolatables, t, x,
+ dof_table);
return extrapolator.getNodalValues();
};
- auto const eval_residuals = [&extrapolator, &local_assemblers,
+ auto const eval_residuals = [num_components, &extrapolator,
+ &local_assemblers,
integration_point_values_method](
- GlobalVector const& /*x*/,
- NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ const double t,
+ 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);
+ extrapolator.calculateResiduals(num_components, extrapolatables, t, x,
+ dof_table);
return extrapolator.getElementResiduals();
};
- return {eval_field, eval_residuals};
+ return {num_components, eval_field, eval_residuals};
}
} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
index 01c7f48fa3612ff60019bcde5e484f6b970cd6ef..b865909f8fda47f575476a7d390f938a62f7aa5f 100644
--- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
+++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
@@ -23,39 +23,75 @@ struct SmallDeformationLocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
virtual std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
};
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index 1ddad16f816b2e012217084111436422be9bf9d8..b89809c8e8e2f9137c744e0947056bb608410c87 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -241,37 +241,55 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 1);
}
std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 3);
}
- std::vector<double> const& getIntPtSigmaYZ(
+ std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
return getIntPtSigma(cache, 4);
}
- std::vector<double> const& getIntPtSigmaXZ(
+ std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -279,30 +297,45 @@ public:
}
std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 0);
}
std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 1);
}
std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 2);
}
std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 3);
}
std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -310,6 +343,9 @@ public:
}
std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
index a1212b69a419f865ac19ef0c47c70e9a0f54454d..df0b44c04b808c5a2d63c7bd6030a8ffd0656adc 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
@@ -69,79 +69,79 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
_nodal_forces->resize(DisplacementDim * mesh.getNumberOfNodes());
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xx", 1,
+ "sigma_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaXX));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yy", 1,
+ "sigma_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaYY));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_zz", 1,
+ "sigma_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaZZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xy", 1,
+ "sigma_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xz", 1,
+ "sigma_xz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaXZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yz", 1,
+ "sigma_yz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtSigmaYZ));
}
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xx", 1,
+ "epsilon_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXX));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yy", 1,
+ "epsilon_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonYY));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_zz", 1,
+ "epsilon_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonZZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xy", 1,
+ "epsilon_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yz", 1,
+ "epsilon_yz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonYZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xz", 1,
+ "epsilon_xz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&SmallDeformationLocalAssemblerInterface::getIntPtEpsilonXZ));
}
}
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 72cfe4c725e63534474bcfbfd378f59d9e16a9b2..ef83e3445ed699e9ecdf88b71a53ac6bccb88378 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -10,6 +10,7 @@
#include "MaterialLib/Adsorption/Adsorption.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "NumLib/Function/Interpolation.h"
@@ -112,7 +113,8 @@ TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
bool is_axially_symmetric,
unsigned const integration_order,
AssemblyParams const& asm_params)
- : _integration_method(integration_order),
+ : _element(e),
+ _integration_method(integration_order),
_shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
IntegrationMethod_, GlobalDim>(
e, is_axially_symmetric, _integration_method)),
@@ -184,18 +186,24 @@ void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::assemble(
template <typename ShapeFunction_, typename IntegrationMethod_,
unsigned GlobalDim>
-std::vector<double> const& TESLocalAssembler<
- ShapeFunction_, IntegrationMethod_,
- GlobalDim>::getIntPtSolidDensity(std::vector<double>& /*cache*/) const
+std::vector<double> const&
+TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+ getIntPtSolidDensity(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/) const
{
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
+std::vector<double> const&
+TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+ getIntPtLoading(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const
{
auto const rho_SR = _d.getData().solid_density;
auto const rho_SR_dry = _d.getAssemblyParameters().rho_SR_dry;
@@ -204,7 +212,7 @@ std::vector<double> const& TESLocalAssembler<
cache.reserve(rho_SR.size());
for (auto const rho : rho_SR) {
- cache.push_back(rho/rho_SR_dry - 1.0);
+ cache.push_back(rho / rho_SR_dry - 1.0);
}
return cache;
@@ -214,7 +222,10 @@ template <typename ShapeFunction_, typename IntegrationMethod_,
unsigned GlobalDim>
std::vector<double> const&
TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
- getIntPtReactionDampingFactor(std::vector<double>& cache) const
+ getIntPtReactionDampingFactor(
+ const double /*t*/, GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& cache) const
{
auto const fac = _d.getData().reaction_adaptor->getReactionDampingFactor();
auto const num_integration_points = _d.getData().solid_density.size();
@@ -227,40 +238,55 @@ TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
template <typename ShapeFunction_, typename IntegrationMethod_,
unsigned GlobalDim>
-std::vector<double> const& TESLocalAssembler<
- ShapeFunction_, IntegrationMethod_,
- GlobalDim>::getIntPtReactionRate(std::vector<double>& /*cache*/) const
+std::vector<double> const&
+TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+ getIntPtReactionRate(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ 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
+std::vector<double> const&
+TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+ getIntPtDarcyVelocity(const double /*t*/,
+ GlobalVector const& current_solution,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<double>& cache) const
{
- return _d.getData().velocity[0];
-}
+ auto const n_integration_points = _integration_method.getNumberOfPoints();
-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];
-}
+ auto const indices = NumLib::getIndices(_element.getID(), dof_table);
+ assert(!indices.empty());
+ auto const local_x = current_solution.get(indices);
+ // local_x is ordered by component, local_x_mat is row major
+ auto const local_x_mat = MathLib::toMatrix<
+ Eigen::Matrix<double, NODAL_DOF, Eigen::Dynamic, Eigen::RowMajor>>(
+ local_x, NODAL_DOF, ShapeFunction_::NPOINTS);
-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];
+ cache.clear();
+ auto cache_mat = MathLib::createZeroedMatrix<
+ Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
+ cache, GlobalDim, n_integration_points);
+
+ for (unsigned i = 0; i < n_integration_points; ++i)
+ {
+ double p, T, x;
+ NumLib::shapeFunctionInterpolate(local_x, _shape_matrices[i].N, p, T,
+ x);
+ const double eta_GR = fluid_viscosity(p, T, x);
+
+ auto const& k = _d.getAssemblyParameters().solid_perm_tensor;
+ cache_mat.col(i).noalias() =
+ k * (_shape_matrices[i].dNdx *
+ local_x_mat.row(COMPONENT_ID_PRESSURE).transpose()) /
+ -eta_GR;
+ }
+
+ return cache;
}
template <typename ShapeFunction_, typename IntegrationMethod_,
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index a69e42750aa9651b568e17da8134255ef718264f..2f47adb3e54a8cee2af01eaddf711a3e1eed28e5 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -32,24 +32,33 @@ public:
std::vector<double> const& local_x_prev_ts) = 0;
virtual std::vector<double> const& getIntPtSolidDensity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtLoading(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtReactionDampingFactor(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtReactionRate(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
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(
+ virtual std::vector<double> const& getIntPtDarcyVelocity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -87,26 +96,38 @@ public:
std::vector<double> const& local_x_prev_ts) override;
std::vector<double> const& getIntPtSolidDensity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override;
std::vector<double> const& getIntPtLoading(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override;
std::vector<double> const& getIntPtReactionDampingFactor(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override;
std::vector<double> const& getIntPtReactionRate(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override;
- std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double> const& getIntPtDarcyVelocity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
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:
+ MeshLib::Element const& _element;
+
IntegrationMethod_ const _integration_method;
std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index db778b975b82bc62ee2f695cc836e2a551a9153e..91ca1a01f375e994a8aa95e650d466cc1967daaf 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -152,17 +152,21 @@ void TESProcess::initializeConcreteProcess(
void TESProcess::initializeSecondaryVariables()
{
// adds a secondary variables to the collection of all secondary variables.
- auto add2nd = [&](std::string const& var_name, unsigned const n_components,
+ auto add2nd = [&](std::string const& var_name,
SecondaryVariableFunctions&& fcts) {
- _secondary_variables.addSecondaryVariable(var_name, n_components,
- std::move(fcts));
+ _secondary_variables.addSecondaryVariable(var_name, std::move(fcts));
};
// creates an extrapolator
- auto makeEx =
- [&](std::vector<double> const& (TESLocalAssemblerInterface::*method)(
- std::vector<double>&)const) -> SecondaryVariableFunctions {
- return ProcessLib::makeExtrapolator(getExtrapolator(),
+ auto makeEx = [&](
+ unsigned const n_components,
+ std::vector<double> const& (TESLocalAssemblerInterface::*method)(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::vector<double>& /*cache*/)
+ const) -> SecondaryVariableFunctions {
+ return ProcessLib::makeExtrapolator(n_components, getExtrapolator(),
_local_assemblers, method);
};
@@ -192,32 +196,28 @@ void TESProcess::initializeSecondaryVariables()
for (auto&& fct : solid_density->getNamedFunctions())
_named_function_caller.addNamedFunction(std::move(fct));
- add2nd("solid_density", 1, solid_density->getExtrapolator());
+ add2nd("solid_density", solid_density->getExtrapolator());
_cached_secondary_variables.emplace_back(std::move(solid_density));
// /////////////////////////////////////////////////////////////////////////
- add2nd("reaction_rate", 1,
- makeEx(&TESLocalAssemblerInterface::getIntPtReactionRate));
+ add2nd("reaction_rate",
+ makeEx(1, &TESLocalAssemblerInterface::getIntPtReactionRate));
- add2nd("velocity_x", 1,
- makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityX));
- if (_mesh.getDimension() >= 2)
- add2nd("velocity_y", 1,
- makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityY));
- if (_mesh.getDimension() >= 3)
- add2nd("velocity_z", 1,
- makeEx(&TESLocalAssemblerInterface::getIntPtDarcyVelocityZ));
+ add2nd("darcy_velocity",
+ makeEx(_mesh.getDimension(),
+ &TESLocalAssemblerInterface::getIntPtDarcyVelocity));
- add2nd("loading", 1, makeEx(&TESLocalAssemblerInterface::getIntPtLoading));
- add2nd("reaction_damping_factor", 1,
- makeEx(&TESLocalAssemblerInterface::getIntPtReactionDampingFactor));
+ add2nd("loading", makeEx(1, &TESLocalAssemblerInterface::getIntPtLoading));
+ add2nd(
+ "reaction_damping_factor",
+ makeEx(1, &TESLocalAssemblerInterface::getIntPtReactionDampingFactor));
- add2nd("relative_humidity", 1,
- {BaseLib::easyBind(&TESProcess::computeRelativeHumidity, this),
+ add2nd("relative_humidity",
+ {1, BaseLib::easyBind(&TESProcess::computeRelativeHumidity, this),
nullptr});
- add2nd("equilibrium_loading", 1,
- {BaseLib::easyBind(&TESProcess::computeEquilibriumLoading, this),
+ add2nd("equilibrium_loading",
+ {1, BaseLib::easyBind(&TESProcess::computeEquilibriumLoading, this),
nullptr});
}
@@ -316,8 +316,8 @@ NumLib::IterationResult TESProcess::postIterationConcreteProcess(
return NumLib::IterationResult::SUCCESS;
}
-GlobalVector const&
-TESProcess::computeVapourPartialPressure(
+GlobalVector const& TESProcess::computeVapourPartialPressure(
+ const double /*t*/,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache)
@@ -348,8 +348,8 @@ TESProcess::computeVapourPartialPressure(
return *result_cache;
}
-GlobalVector const&
-TESProcess::computeRelativeHumidity(
+GlobalVector const& TESProcess::computeRelativeHumidity(
+ double const /*t*/,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache)
@@ -385,8 +385,8 @@ TESProcess::computeRelativeHumidity(
return *result_cache;
}
-GlobalVector const&
-TESProcess::computeEquilibriumLoading(
+GlobalVector const& TESProcess::computeEquilibriumLoading(
+ double const /*t*/,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache)
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index f75afb5b181adae3720af453c4eccd5df3f5f1d4..87449a97678dea23577fda6ed7669f9efe20634a 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -69,16 +69,19 @@ private:
StaggeredCouplingTerm const& coupling_term) override;
GlobalVector const& computeVapourPartialPressure(
+ const double t,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache);
GlobalVector const& computeRelativeHumidity(
+ const double t,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache);
GlobalVector const& computeEquilibriumLoading(
+ const double t,
GlobalVector const& x,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::unique_ptr<GlobalVector>& result_cache);
diff --git a/ProcessLib/TES/Tests.cmake b/ProcessLib/TES/Tests.cmake
index bdcbc6ee8e75f8bfcfdbc84ce822190f9d7028c6..57d9924d0294140fc1412742c732602b9526f57d 100644
--- a/ProcessLib/TES/Tests.cmake
+++ b/ProcessLib/TES/Tests.cmake
@@ -12,7 +12,6 @@ AddTest(
tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000.vtu temperature temperature
tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000.vtu vapour_partial_pressure vapour_partial_pressure
tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000.vtu solid_density solid_density
- tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_pcs_0_ts_19_t_0.100000.vtu solid_density fct_solid_density
)
AddTest(
@@ -28,7 +27,6 @@ AddTest(
tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000.vtu temperature temperature
tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000.vtu vapour_partial_pressure vapour_partial_pressure
tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000.vtu solid_density solid_density
- tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000.vtu solid_density fct_solid_density
)
AddTest(
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
index f9a9d2eb201cc8f5983dad69ec6207d2de449223..fe07bb3dbd587d0993df4a410bbb82e27fc1fed3 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
@@ -59,9 +59,15 @@ class ThermalTwoPhaseFlowWithPPLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtWettingPressure(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -138,6 +144,9 @@ public:
}
std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_saturation.empty());
@@ -145,6 +154,9 @@ public:
}
std::vector<double> const& getIntPtWettingPressure(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_pressure_wetting.empty());
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
index 25733bc90bd933d94f5dcd9ad6e796c0e193c05c..8cf9ac87fde8b53d50c809e1b4d1370f2d3a6f77 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPProcess.cpp
@@ -57,14 +57,14 @@ void ThermalTwoPhaseFlowWithPPProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "saturation", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "saturation",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ThermalTwoPhaseFlowWithPPLocalAssemblerInterface::
getIntPtSaturation));
_secondary_variables.addSecondaryVariable(
- "pressure_wetting", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "pressure_wetting",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ThermalTwoPhaseFlowWithPPLocalAssemblerInterface::
getIntPtWettingPressure));
}
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index 1ee2f46107ec206c6eabec720428f3d9904f5d82..dd85ee6c7a3bd8c1f99337d7025324018d90ed3a 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -102,39 +102,75 @@ struct ThermoMechanicsLocalAssemblerInterface
public NumLib::ExtrapolatableElement
{
virtual std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
};
@@ -385,30 +421,45 @@ public:
}
std::vector<double> const& getIntPtSigmaXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 0);
}
std::vector<double> const& getIntPtSigmaYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 1);
}
std::vector<double> const& getIntPtSigmaZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 2);
}
std::vector<double> const& getIntPtSigmaXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtSigma(cache, 3);
}
std::vector<double> const& getIntPtSigmaYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -416,6 +467,9 @@ public:
}
std::vector<double> const& getIntPtSigmaXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -423,30 +477,45 @@ public:
}
std::vector<double> const& getIntPtEpsilonXX(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 0);
}
std::vector<double> const& getIntPtEpsilonYY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 1);
}
std::vector<double> const& getIntPtEpsilonZZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 2);
}
std::vector<double> const& getIntPtEpsilonXY(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
return getIntPtEpsilon(cache, 3);
}
std::vector<double> const& getIntPtEpsilonYZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
@@ -454,6 +523,9 @@ public:
}
std::vector<double> const& getIntPtEpsilonXZ(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
assert(DisplacementDim == 3);
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
index 7e6c4c290da5fbc92f3f943623cde64425925e9f..abf3837395987b822a3b4d29bd01602b80bfec83 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsProcess.h
@@ -79,77 +79,77 @@ private:
NumLib::ComponentOrder::BY_LOCATION));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xx", 1,
+ "sigma_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXX));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yy", 1,
+ "sigma_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYY));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_zz", 1,
+ "sigma_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaZZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xy", 1,
+ "sigma_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "sigma_xz", 1,
+ "sigma_xz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaXZ));
Base::_secondary_variables.addSecondaryVariable(
- "sigma_yz", 1,
+ "sigma_yz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtSigmaYZ));
}
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xx", 1,
+ "epsilon_xx",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXX));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yy", 1,
+ "epsilon_yy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonYY));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_zz", 1,
+ "epsilon_zz",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonZZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xy", 1,
+ "epsilon_xy",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::getIntPtEpsilonXY));
if (DisplacementDim == 3)
{
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_yz", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "epsilon_yz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::
getIntPtEpsilonYZ));
Base::_secondary_variables.addSecondaryVariable(
- "epsilon_xz", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "epsilon_xz",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&ThermoMechanicsLocalAssemblerInterface::
getIntPtEpsilonXZ));
}
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
index 9c53d8678c665a555f84731e16705ff95c8e4db1..e532c118113501d0573ae22c026d445bd8873084 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
@@ -59,9 +59,15 @@ class TwoPhaseFlowWithPPLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtWetPressure(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -135,6 +141,9 @@ public:
}
std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_saturation.empty());
@@ -142,6 +151,9 @@ public:
}
std::vector<double> const& getIntPtWetPressure(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_pressure_wet.empty());
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 6c855e3ba4cb4ce6e64cea3ecbdfb6ad9e519736..6e7fb5f328c461919811627d356fc4d4822b770a 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -56,16 +56,16 @@ void TwoPhaseFlowWithPPProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "saturation", 1,
+ "saturation",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&TwoPhaseFlowWithPPLocalAssemblerInterface::getIntPtSaturation));
_secondary_variables.addSecondaryVariable(
- "pressure_wet", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
- &TwoPhaseFlowWithPPLocalAssemblerInterface::
- getIntPtWetPressure));
+ "pressure_wet",
+ makeExtrapolator(
+ 1, getExtrapolator(), _local_assemblers,
+ &TwoPhaseFlowWithPPLocalAssemblerInterface::getIntPtWetPressure));
}
void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
index 28c994d8f84c7f46d78b9e159c3f7f802d855aa0..b505b7721eda9a858bfa52413fa343f5bf85156c 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
@@ -61,9 +61,15 @@ class TwoPhaseFlowWithPrhoLocalAssemblerInterface
{
public:
virtual std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtNonWettingPressure(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
};
@@ -138,6 +144,9 @@ public:
}
std::vector<double> const& getIntPtSaturation(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_saturation.empty());
@@ -145,6 +154,9 @@ public:
}
std::vector<double> const& getIntPtNonWettingPressure(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(!_pressure_nonwetting.empty());
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
index e57594d1c2035a383ae824e116fd4fd8226ac3d0..ef8105649ef78cb1d8637f5d45a7d8aa5b720c40 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -56,14 +56,14 @@ void TwoPhaseFlowWithPrhoProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable(
- "saturation", 1,
+ "saturation",
makeExtrapolator(
- getExtrapolator(), _local_assemblers,
+ 1, getExtrapolator(), _local_assemblers,
&TwoPhaseFlowWithPrhoLocalAssemblerInterface::getIntPtSaturation));
_secondary_variables.addSecondaryVariable(
- "pressure_nonwetting", 1,
- makeExtrapolator(getExtrapolator(), _local_assemblers,
+ "pressure_nonwetting",
+ makeExtrapolator(1, getExtrapolator(), _local_assemblers,
&TwoPhaseFlowWithPrhoLocalAssemblerInterface::
getIntPtNonWettingPressure));
}
diff --git a/Tests/Data b/Tests/Data
index 90cf6764663bf30b32e89e72737cb3413d081609..92a819ff3d83ffd9d77517c071a756af88ebbd47 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 90cf6764663bf30b32e89e72737cb3413d081609
+Subproject commit 92a819ff3d83ffd9d77517c071a756af88ebbd47
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index a17d23fa0e2c367ffe1d789093e764531ba8eb3d..075556d43ce5aec93cdad88f370e124ccfd2264e 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -55,14 +55,25 @@ public:
std::vector<double> const& local_nodal_values) = 0;
virtual std::vector<double> const& getStoredQuantity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getDerivedQuantity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const = 0;
};
using IntegrationPointValuesMethod = std::vector<double> const& (
- LocalAssemblerDataInterface::*)(std::vector<double>&)const;
+ LocalAssemblerDataInterface::*)(const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::
+ LocalToGlobalIndexMap const& /*dof_table*/
+ ,
+ std::vector<double>& /*cache*/) const;
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
@@ -95,12 +106,18 @@ public:
}
std::vector<double> const& getStoredQuantity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
return _int_pt_values;
}
std::vector<double> const& getDerivedQuantity(
+ const double /*t*/,
+ GlobalVector const& /*current_solution*/,
+ NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& cache) const override
{
cache.clear();
@@ -169,13 +186,15 @@ public:
}
std::pair<GlobalVector const*, GlobalVector const*> extrapolate(
- IntegrationPointValuesMethod method) const
+ IntegrationPointValuesMethod method, const double t,
+ const GlobalVector& x) const
{
auto const extrapolatables =
NumLib::makeExtrapolatable(_local_assemblers, method);
- _extrapolator->extrapolate(extrapolatables);
- _extrapolator->calculateResiduals(extrapolatables);
+ _extrapolator->extrapolate(1, extrapolatables, t, x, *_dof_table);
+ _extrapolator->calculateResiduals(1, extrapolatables, t, x,
+ *_dof_table);
return {&_extrapolator->getNodalValues(),
&_extrapolator->getElementResiduals()};
@@ -201,7 +220,10 @@ void extrapolate(TestProcess const& pcs, IntegrationPointValuesMethod method,
auto const tolerance_dx = 30.0 * std::numeric_limits<double>::epsilon();
auto const tolerance_res = 15.0 * std::numeric_limits<double>::epsilon();
- auto const result = pcs.extrapolate(method);
+ const double t = 0.0;
+
+ auto const result =
+ pcs.extrapolate(method, t, expected_extrapolated_global_nodal_values);
auto const& x_extra = *result.first;
auto const& residual = *result.second;