ufz · chleh · Aug 17, 2016 · Jun 23, 2016 · Jul 7, 2016 · Jul 7, 2016
diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
@@ -198,6 +198,7 @@ if(NOT OGS_USE_MPI)
         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(
@@ -213,6 +214,7 @@ if(NOT OGS_USE_MPI)
         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
     )
 
 else()

diff --git a/MathLib/LinAlg/MatrixVectorTraits.h b/MathLib/LinAlg/MatrixVectorTraits.h
@@ -10,14 +10,13 @@
 #ifndef MATHLIB_MATRIX_VECTOR_TRAITS_H
 #define MATHLIB_MATRIX_VECTOR_TRAITS_H
 
-#include<memory>
+#include <memory>
+#include "MatrixSpecifications.h"
 
 namespace MathLib
 {
 template<typename Matrix>
 struct MatrixVectorTraits;
-
-struct MatrixSpecifications;
 }
 
 #define SPECIALIZE_MATRIX_VECTOR_TRAITS(MATVEC, IDX) \

diff --git a/NumLib/Extrapolation/ExtrapolatableElementCollection.h b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
@@ -10,6 +10,7 @@
 #ifndef NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENTCOLLECTION_H
 #define NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENTCOLLECTION_H
 
+#include <vector>
 #include "ExtrapolatableElement.h"
 
 namespace NumLib

diff --git a/ProcessLib/CachedSecondaryVariable.cpp b/ProcessLib/CachedSecondaryVariable.cpp
@@ -0,0 +1,70 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CachedSecondaryVariable.h"
+#include "BaseLib/Functional.h"
+#include "MathLib/LinAlg/LinAlg.h"
+
+namespace ProcessLib
+{
+std::vector<NumLib::NamedFunction> CachedSecondaryVariable::getNamedFunctions()
+    const
+{
+    return {{_internal_variable_name,
+             std::vector<std::string>{},
+             BaseLib::easyBind(&CachedSecondaryVariable::getValue, this)}};
+}
+
+double CachedSecondaryVariable::getValue() const
+{
+    if (_needs_recomputation)
+        evalFieldNoArgs();
+    return _cached_nodal_values.get(_context.index);
+}
+
+SecondaryVariableFunctions CachedSecondaryVariable::getExtrapolator()
+{
+    // TODO copied from makeExtrapolator()
+    auto const eval_residuals = [this](
+        GlobalVector const& /*x*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::unique_ptr<GlobalVector> & /*result_cache*/
+        ) -> GlobalVector const& {
+        _extrapolator.calculateResiduals(*_extrapolatables);
+        return _extrapolator.getElementResiduals();
+    };
+    return {BaseLib::easyBind(&CachedSecondaryVariable::evalField, this),
+            eval_residuals};
+}
+
+GlobalVector const& CachedSecondaryVariable::evalField(
+    GlobalVector const& /*x*/,
+    NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+    std::unique_ptr<GlobalVector>& /*result_cache*/
+    ) const
+{
+    return evalFieldNoArgs();
+}
+
+GlobalVector const& CachedSecondaryVariable::evalFieldNoArgs() const
+{
+    if (!_needs_recomputation) {
+        DBUG("%s does not need to be recomputed. Returning cached values",
+             _internal_variable_name.c_str());
+        return _cached_nodal_values;
+    }
+    DBUG("Recomputing %s.", _internal_variable_name.c_str());
+    _extrapolator.extrapolate(*_extrapolatables);
+    auto const& nodal_values = _extrapolator.getNodalValues();
+    MathLib::LinAlg::copy(nodal_values, _cached_nodal_values);
+    _needs_recomputation = false;
+    return nodal_values;
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/CachedSecondaryVariable.h b/ProcessLib/CachedSecondaryVariable.h
@@ -0,0 +1,91 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_CACHEDSECONDARYVARIABLE_H
+#define PROCESSLIB_CACHEDSECONDARYVARIABLE_H
+
+#include "NumLib/Extrapolation/ExtrapolatableElementCollection.h"
+#include "NumLib/NamedFunctionProvider.h"
+#include "NumLib/NumericsConfig.h"
+#include "SecondaryVariable.h"
+#include "SecondaryVariableContext.h"
+
+namespace ProcessLib
+{
+/*! Secondary variable which is extrapolated from integration points to mesh
+ * nodes; the resulting extrapolated values are cached for subsequent use.
+ */
+class CachedSecondaryVariable final : public NumLi
1E0A
b::NamedFunctionProvider
+{
+public:
+    /*! Constructs a new instance.
+     *
+     * \param internal_variable_name the variable's name
+     * \param extrapolator extrapolates integration point values to nodal
+     * values.
+     * \param local_assemblers provide the integration point values
+     * \param integration_point_values_method extracts the integration point
+     * values from the \c local_assemblers
+     * \param context needed s.t. this class can act as a NamedFunction
+     */
+    template <typename LocalAssemblerCollection,
              typename IntegrationPointValuesMethod>
+    CachedSecondaryVariable(
+        std::string const& internal_variable_name,
+        NumLib::Extrapolator& extrapolator,
+        LocalAssemblerCollection const& local_assemblers,
+        IntegrationPointValuesMethod integration_point_values_method,
+        SecondaryVariableContext const& context)
+        : _extrapolator(extrapolator)
+        , _extrapolatables(new NumLib::ExtrapolatableLocalAssemblerCollection<
+                           LocalAssemblerCollection>{
+              local_assemblers, integration_point_values_method})
+        , _context(context)
+        , _internal_variable_name(internal_variable_name)
+    {
+    }
+
+    CachedSecondaryVariable(CachedSecondaryVariable const&) = delete;
+    CachedSecondaryVariable(CachedSecondaryVariable&&) = delete;
+
+    std::vector<NumLib::NamedFunction> getNamedFunctions() const override;
+
+    //! Returns extrapolation functions that compute the secondary variable.
+    SecondaryVariableFunctions getExtrapolator();
+
+    //! Set that recomputation is necessary.
+    void expire() { _needs_recomputation = true; }
+
+private:
+    //! Provides the value at the current index of the _context.
+    double getValue() const;
+
+    //! Computes the secondary Variable.
+    GlobalVector const& evalField(
+        GlobalVector const& /*x*/,
+        NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+        std::unique_ptr<GlobalVector>& /*result_cache*/
+        ) const;
+
+    //! Computes the secondary Variable.
+    GlobalVector const& evalFieldNoArgs() const;
+
+    //! Cache for the computed values.
+    mutable GlobalVector _cached_nodal_values;
+    mutable bool _needs_recomputation = true;
+
+    NumLib::Extrapolator& _extrapolator;
+    std::unique_ptr<NumLib::ExtrapolatableElementCollection> _extrapolatables;
+    SecondaryVariableContext const& _context;
+    std::string const _internal_variable_name;
+};
+
+}  // namespace ProcessLib
+
+#endif // PROCESSLIB_CACHEDSECONDARYVARIABLE_H
diff --git a/ProcessLib/GlobalVectorFromNamedFunction.cpp b/ProcessLib/GlobalVectorFromNamedFunction.cpp
@@ -0,0 +1,63 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "GlobalVectorFromNamedFunction.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "ProcessLib/Utils/VectorUtil.h"
+
+namespace ProcessLib
+{
+GlobalVectorFromNamedFunction::GlobalVectorFromNamedFunction(
+    NumLib::SpecificFunctionCaller&& function_caller,
+    MeshLib::Mesh const& mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table_single,
+    SecondaryVariableContext& context)
+    : _function_caller(std::move(function_caller))
+    , _mesh(mesh)
+    , _dof_table_single(dof_table_single)
+    , _context(context)
+{
+    assert(dof_table_single.getNumberOfComponents() == 1);
+}
+
+GlobalVector const& GlobalVectorFromNamedFunction::call(
+    GlobalVector const& x,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::unique_ptr<GlobalVector>& result)
+{
+    result = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+        {_dof_table_single.dofSizeWithoutGhosts(),
+         _dof_table_single.dofSizeWithoutGhosts(),
+         &_dof_table_single.getGhostIndices(), nullptr});
+
+    GlobalIndexType nnodes = _mesh.getNumberOfNodes();
+
+    auto const n_args = _function_caller.getNumberOfUnboundArguments();
+    assert(dof_table.getNumberOfComponents() == n_args);
+    std::vector<double> args(n_args);
+
+    for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id) {
+        // TODO maybe fill args via callback mechanism or remove this class entirely.
+        // Caution: The order of args will be the same as the order of the
+        // components in the global vector!
+        for (std::size_t i = 0; i < n_args; ++i) {
+            args[i] = getNodalValue(x, _mesh, dof_table, node_id, i);
+        }
+
+        _context.index = node_id;
+        auto const value = _function_caller.call(args);
+
+        // TODO Problems with PETSc? (local vs. global index)
+        result->set(node_id, value);
+    }
+
+    return *result;
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/GlobalVectorFromNamedFunction.h b/ProcessLib/GlobalVectorFromNamedFunction.h
@@ -0,0 +1,58 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef PROCESSLIB_GLOBALVECTORFROMNAMEDFUNCTION_H
+#define PROCESSLIB_GLOBALVECTORFROMNAMEDFUNCTION_H
+
+#include "MeshLib/Mesh.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/NumericsConfig.h"
+#include "NumLib/NamedFunctionCaller.h"
+#include "SecondaryVariableContext.h"
+
+namespace ProcessLib
+{
+//! Computes a global vector from a NamedFunction (which can only compute double
+//! values).
+class GlobalVectorFromNamedFunction final
+{
+public:
+    /*! Constructs a new instance.
+     *
+     * \param function_caller will provide the individual entries of the
+     * GlobalVector to be computed.
+     * \param mesh to which the \c dof_table_single is associated
+     * \param dof_table_single used for constructing the GlobalVector
+     * \param context used by the \c function_caller to access "auxiliary
+     * unbound arguments".
+     */
+    GlobalVectorFromNamedFunction(
+        NumLib::SpecificFunctionCaller&&
10000
amp; function_caller,
+        MeshLib::Mesh const& mesh,
+        NumLib::LocalToGlobalIndexMap const& dof_table_single,
+        SecondaryVariableContext& context);
+
+    //! Computes the GlobalVector.
+    //!
+    //! 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,
+                             NumLib::LocalToGlobalIndexMap const& dof_table,
+                             std::unique_ptr<GlobalVector>& result);
+
+private:
+    NumLib::SpecificFunctionCaller _function_caller;
+    MeshLib::Mesh const& _mesh;
+    NumLib::LocalToGlobalIndexMap const& _dof_table_single;
+    SecondaryVariableContext& _context;
+};
+}  // namespace ProcessLib
+
+#endif  // PROCESSLIB_GLOBALVECTORFROMNAMEDFUNCTION_H
diff --git a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp
@@ -9,6 +9,7 @@
 
 #include "CreateGroundwaterFlowProcess.h"
 
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
 #include "GroundwaterFlowProcess.h"
 #include "GroundwaterFlowProcessData.h"
 
@@ -48,25 +49,22 @@ std::unique_ptr<Process> createGroundwaterFlowProcess(
 
     GroundwaterFlowProcessData process_data{hydraulic_conductivity};
 
-    SecondaryVariableCollection secondary_variables{
-        //! \ogs_file_param{process__secondary_variables}
-        config.getConfigSubtreeOptional("secondary_variables"),
-        {//! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_x}
-         "darcy_velocity_x",
-         //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_y}
-         "darcy_velocity_y",
-         //! \ogs_file_param_special{process__GROUNDWATER_FLOW__secondary_variables__darcy_velocity_z}
-         "darcy_velocity_z"}};
+    SecondaryVariableCollection secondary_variables;
 
-    ProcessOutput
-        //! \ogs_file_param{process__output}
-        process_output{config.getConfigSubtree("output"), process_variables,
-                       secondary_variables};
+    NumLib::NamedFunctionCaller named_function_caller(
+        {"GWFlow_pressure"});
+
+    ProcessLib::parseSecondaryVariables(config, secondary_variables,
+                                        named_function_caller);
+
+    //! \ogs_file_param{process__output}
+    ProcessOutput process_output{config.getConfigSubtree("output")};
 
     return std::unique_ptr<Process>{new GroundwaterFlowProcess{
         mesh, nonlinear_solver, std::move(time_discretization),
         std::move(process_variables), std::move(process_data),
-        std::move(secondary_variables), std::move(process_output)}};
+        std::move(secondary_variables), std::move(process_output),
+        std::move(named_function_caller)}};
 }
 
 }  // namespace GroundwaterFlow

diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -24,10 +24,11 @@ GroundwaterFlowProcess::GroundwaterFlowProcess(
     std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
     GroundwaterFlowProcessData&& process_data,
     SecondaryVariableCollection&& secondary_variables,
-    ProcessOutput&& process_output)
+    ProcessOutput&& process_output,
+    NumLib::NamedFunctionCaller&& named_function_caller)
     : Process(mesh, nonlinear_solver, std::move(time_discretization),
               std::move(process_variables), std::move(secondary_variables),
-              std::move(process_output)),
+              std::move(process_output), std::move(named_function_caller)),
       _process_data(std::move(process_data))
 {
     if (dynamic_cast<NumLib::ForwardEuler*>(

diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -33,7 +33,8 @@ class GroundwaterFlowProcess final : public Process
             process_variables,
         GroundwaterFlowProcessData&& process_data,
         SecondaryVariableCollection&& secondary_variables,
-        ProcessOutput&& process_output);
+        ProcessOutput&& process_output,
+        NumLib::NamedFunctionCaller&& named_function_caller);
 
     //! \name ODESystem interface
     //! @{