diff --git a/BaseLib/uniqueInsert.h b/BaseLib/uniqueInsert.h
index 5bbdbf7fe37d06976dee0defe94b20a7abbeaeb0..047509e9e7250ec9d6a829e9eea7b5e552d46da7 100644
--- a/BaseLib/uniqueInsert.h
+++ b/BaseLib/uniqueInsert.h
@@ -44,6 +44,31 @@ void insertIfKeyUniqueElseError(
}
}
+//! Inserts the given \c key with the given \c value into the \c map if neither an entry
+//! with the given \c key nor an entry with the given \c value already exists;
+//! otherwise an \c error_message is printed and the program is aborted.
+template<typename Map, typename Key, typename Value>
+void insertIfKeyValueUniqueElseError(
+ Map& map, Key const& key, Value&& value,
+ std::string const& error_message)
+{
+ auto value_compare = [&value](typename Map::value_type const& elem) {
+ return value == elem.second;
+ };
+
+ if (std::find_if(map.cbegin(), map.cend(), value_compare) != map.cend())
+ {
+ ERR("%s Value `%s' already exists.", error_message.c_str(), tostring(value).c_str());
+ std::abort();
+ }
+
+ auto const inserted = map.emplace(key, std::forward<Value>(value));
+ if (!inserted.second) { // insertion failed, i.e., key already exists
+ ERR("%s Key `%s' already exists.", error_message.c_str(), tostring(key).c_str());
+ std::abort();
+ }
+}
+
//! Returns the value of \c key from the given \c map if such an entry exists;
//! otherwise an \c error_message is printed and the program is aborted.
//! Cf. also the const overload below.
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index d686c625dce35d738e9cdcadb17c7dce9ab3347d..b177164304962a752a6c66014b7570c7c60a295c 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
@@ -52,7 +52,11 @@ public:
MathLib::MatrixSpecifications const& matrix_specs)
: _nodal_values(MathLib::GlobalVectorProvider<GlobalVector>::provider.
getVector(matrix_specs))
+#ifndef USE_PETSC
, _residuals(matrix_specs.dof_table->size())
+#else
+ , _residuals(matrix_specs.dof_table->size(), false)
+#endif
, _local_to_global(*matrix_specs.dof_table)
{}
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index d41e69d5ca79c55b91244afa463a0be0dc0917a5..331bad5a960439e03dc550aaaad0f2ccf4e9caf6 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -26,14 +26,27 @@ namespace ProcessLib
namespace GroundwaterFlow
{
+enum class IntegrationPointValue {
+ DarcyVelocityX,
+ DarcyVelocityY,
+ DarcyVelocityZ
+};
+
const unsigned NUM_NODAL_DOF = 1;
+template <typename GlobalMatrix, typename GlobalVector>
+class GroundwaterFlowLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface<GlobalMatrix, GlobalVector>
+ , public NumLib::Extrapolatable<GlobalVector, IntegrationPointValue>
+{};
+
template <typename ShapeFunction,
typename IntegrationMethod,
typename GlobalMatrix,
typename GlobalVector,
unsigned GlobalDim>
-class LocalAssemblerData : public ProcessLib::LocalAssemblerInterface<GlobalMatrix, GlobalVector>
+class LocalAssemblerData
+ : public GroundwaterFlowLocalAssemblerInterface<GlobalMatrix, GlobalVector>
{
using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
@@ -64,7 +77,7 @@ public:
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
}
- void assemble(double const /*t*/, std::vector<double> const& /*local_x*/) override
+ void assemble(double const /*t*/, std::vector<double> const& local_x) override
{
_localA.setZero();
_localRhs.setZero();
@@ -72,13 +85,23 @@ public:
IntegrationMethod integration_method(_integration_order);
unsigned const n_integration_points = integration_method.getNPoints();
- for (std::size_t ip(0); ip < n_integration_points; ip++) {
+ for (std::size_t ip(0); ip < n_integration_points; ip++)
+ {
auto const& sm = _shape_matrices[ip];
auto const& wp = integration_method.getWeightedPoint(ip);
-
auto const k = _process_data.hydraulic_conductivity(_element);
+
_localA.noalias() += sm.dNdx.transpose() * k * sm.dNdx *
sm.detJ * wp.getWeight();
+
+ // Darcy velocity only computed for output.
+ auto const darcy_velocity = (k * sm.dNdx *
+ Eigen::Map<const NodalVectorType>(local_x.data(), ShapeFunction::NPOINTS)
+ ).eval();
+
+ for (unsigned d=0; d<GlobalDim; ++d) {
+ _darcy_velocities[d][ip] = darcy_velocity[d];
+ }
}
}
@@ -90,6 +113,34 @@ public:
b.add(indices.rows, _localRhs);
}
+ Eigen::Map<const Eigen::VectorXd>
+ getShapeMatrix(const unsigned integration_point) const override
+ {
+ auto const& N = _shape_matrices[integration_point].N;
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::VectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const&
+ getIntegrationPointValues(IntegrationPointValue const property,
+ std::vector<double>& /*cache*/) const override
+ {
+ switch (property)
+ {
+ case IntegrationPointValue::DarcyVelocityX:
+ return _darcy_velocities[0];
+ case IntegrationPointValue::DarcyVelocityY:
+ assert(GlobalDim > 1);
+ return _darcy_velocities[1];
+ case IntegrationPointValue::DarcyVelocityZ:
+ assert(GlobalDim > 2);
+ return _darcy_velocities[2];
+ }
+
+ std::abort();
+ }
+
private:
MeshLib::Element const& _element;
std::vector<ShapeMatrices> _shape_matrices;
@@ -99,6 +150,10 @@ private:
NodalVectorType _localRhs;
unsigned const _integration_order;
+
+ std::vector<std::vector<double>> _darcy_velocities
+ = std::vector<std::vector<double>>(
+ GlobalDim, std::vector<double>(ShapeFunction::NPOINTS));
};
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index 94635b3b10cfc3aaa3e592b8bef7c6591871554d..603c4978c4db0de2bf8d155ec1964a643a1f041c 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -13,20 +13,16 @@
#include <cassert>
#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
#include "ProcessLib/Process.h"
#include "ProcessLib/Utils/CreateLocalAssemblers.h"
#include "GroundwaterFlowFEM.h"
#include "GroundwaterFlowProcessData.h"
-namespace MeshLib
-{
- class Element;
-}
namespace ProcessLib
{
-
namespace GroundwaterFlow
{
@@ -44,9 +40,14 @@ public:
typename Base::NonlinearSolver& nonlinear_solver,
std::unique_ptr<typename Base::TimeDiscretization>&& time_discretization,
std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
- GroundwaterFlowProcessData&& process_data)
+ GroundwaterFlowProcessData&& process_data,
+ SecondaryVariableCollection<GlobalVector>&& secondary_variables,
+ ProcessOutput<GlobalVector>&& process_output
+ )
: Process<GlobalSetup>(mesh, nonlinear_solver, std::move(time_discretization),
- std::move(process_variables))
+ std::move(process_variables),
+ std::move(secondary_variables),
+ std::move(process_output))
, _process_data(std::move(process_data))
{
if (dynamic_cast<NumLib::ForwardEuler<GlobalVector>*>(
@@ -75,15 +76,21 @@ public:
private:
using LocalAssemblerInterface =
- ProcessLib::LocalAssemblerInterface<GlobalMatrix, GlobalVector>;
+ GroundwaterFlowLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
using GlobalAssembler = AssemblerLib::VectorMatrixAssembler<
GlobalMatrix, GlobalVector, LocalAssemblerInterface,
NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
- void createAssemblers(AssemblerLib::LocalToGlobalIndexMap const& dof_table,
- MeshLib::Mesh const& mesh,
- unsigned const integration_order) override
+ using ExtrapolatorInterface = NumLib::Extrapolator<
+ GlobalVector, IntegrationPointValue, LocalAssemblerInterface>;
+ using ExtrapolatorImplementation = NumLib::LocalLinearLeastSquaresExtrapolator<
+ GlobalVector, IntegrationPointValue, LocalAssemblerInterface>;
+
+ void initializeConcreteProcess(
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order) override
{
DBUG("Create global assembler.");
_global_assembler.reset(new GlobalAssembler(dof_table));
@@ -92,6 +99,28 @@ private:
mesh.getDimension(), mesh.getElements(),
dof_table, integration_order, _local_assemblers,
_process_data);
+
+ // TOOD Later on the DOF table can change during the simulation!
+ _extrapolator.reset(
+ new ExtrapolatorImplementation(Base::getMatrixSpecifications()));
+
+ Base::_secondary_variables.addSecondaryVariable(
+ "darcy_velocity_x", 1,
+ makeExtrapolator(IntegrationPointValue::DarcyVelocityX, *_extrapolator,
+ _local_assemblers));
+
+ if (mesh.getDimension() > 1) {
+ Base::_secondary_variables.addSecondaryVariable(
+ "darcy_velocity_y", 1,
+ makeExtrapolator(IntegrationPointValue::DarcyVelocityY, *_extrapolator,
+ _local_assemblers));
+ }
+ if (mesh.getDimension() > 2) {
+ Base::_secondary_variables.addSecondaryVariable(
+ "darcy_velocity_z", 1,
+ makeExtrapolator(IntegrationPointValue::DarcyVelocityZ, *_extrapolator,
+ _local_assemblers));
+ }
}
void assembleConcreteProcess(const double t, GlobalVector const& x,
@@ -110,6 +139,8 @@ private:
std::unique_ptr<GlobalAssembler> _global_assembler;
std::vector<std::unique_ptr<LocalAssemblerInterface>> _local_assemblers;
+
+ std::unique_ptr<ExtrapolatorInterface> _extrapolator;
};
template <typename GlobalSetup>
@@ -142,11 +173,21 @@ createGroundwaterFlowProcess(
hydraulic_conductivity
};
+ SecondaryVariableCollection<typename GlobalSetup::VectorType>
+ secondary_variables{config.getConfSubtreeOptional("secondary_variables"),
+ { "darcy_velocity_x", "darcy_velocity_y", "darcy_velocity_z" }};
+
+ ProcessOutput<typename GlobalSetup::VectorType>
+ process_output{config.getConfSubtree("output"),
+ process_variables, secondary_variables};
+
return std::unique_ptr<GroundwaterFlowProcess<GlobalSetup>>{
new GroundwaterFlowProcess<GlobalSetup>{
mesh, nonlinear_solver,std::move(time_discretization),
std::move(process_variables),
- std::move(process_data)
+ std::move(process_data),
+ std::move(secondary_variables),
+ std::move(process_output)
}
};
}
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index cdf0fe4f329d4626e35ee5c1f10042f5eb30ae23..8e14d4cc9cda7eb8df156b58da2e73ba824bea6b 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -16,8 +16,18 @@
#include "NumLib/ODESolver/TimeDiscretization.h"
#include "NumLib/ODESolver/NonlinearSolver.h"
+#include "DirichletBc.h"
+#include "NeumannBc.h"
+#include "NeumannBcAssembler.h"
#include "Parameter.h"
+#include "ProcessOutput.h"
#include "ProcessVariable.h"
+#include "UniformDirichletBoundaryCondition.h"
+
+namespace MeshLib
+{
+class Mesh;
+}
namespace ProcessLib
{
@@ -37,15 +47,20 @@ public:
using NonlinearSolver = NumLib::NonlinearSolverBase<GlobalMatrix, GlobalVector>;
using TimeDiscretization = NumLib::TimeDiscretization<GlobalVector>;
- Process(MeshLib::Mesh& mesh,
- NonlinearSolver& nonlinear_solver,
- std::unique_ptr<TimeDiscretization>&& time_discretization,
- std::vector<std::reference_wrapper<ProcessVariable>>&&
- process_variables)
+ Process(
+ MeshLib::Mesh& mesh,
+ NonlinearSolver& nonlinear_solver,
+ std::unique_ptr<TimeDiscretization>&& time_discretization,
+ std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ SecondaryVariableCollection<GlobalVector>&& secondary_variables,
+ ProcessOutput<GlobalVector>&& process_output
+ )
: _mesh(mesh)
, _nonlinear_solver(nonlinear_solver)
, _time_discretization(std::move(time_discretization))
, _process_variables(std::move(process_variables))
+ , _process_output(std::move(process_output))
+ , _secondary_variables(std::move(secondary_variables))
{}
/// Preprocessing before starting assembly for new timestep.
@@ -61,48 +76,8 @@ public:
const unsigned /*timestep*/,
GlobalVector const& x) const
{
- DBUG("Process output.");
-
- // Copy result
-#ifdef USE_PETSC
- // TODO It is also possible directly to copy the data for single process
- // variable to a mesh property. It needs a vector of global indices and
- // some PETSc magic to do so.
- std::vector<double> x_copy(x.getLocalSize() + x.getGhostSize());
-#else
- std::vector<double> x_copy(x.size());
-#endif
- x.copyValues(x_copy);
-
- std::size_t const n_nodes = _mesh.getNNodes();
- for (ProcessVariable& pv : _process_variables)
- {
- auto& output_data = pv.getOrCreateMeshProperty();
-
- int const n_components = pv.getNumberOfComponents();
- for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
- {
- MeshLib::Location const l(_mesh.getID(),
- MeshLib::MeshItemType::Node, node_id);
- // TODO extend component ids to multiple process variables.
- for (int component_id = 0; component_id < n_components;
- ++component_id)
- {
- auto const index =
- _local_to_global_index_map->getLocalIndex(
- l, component_id, x.getRangeBegin(),
- x.getRangeEnd());
-
- output_data[node_id * n_components + component_id] =
- x_copy[index];
- }
- }
- }
-
- // Write output file
- DBUG("Writing output to \'%s\'.", file_name.c_str());
- MeshLib::IO::VtuInterface vtu_interface(&_mesh, vtkXMLWriter::Binary, true);
- vtu_interface.writeToFile(file_name);
+ doProcessOutput(file_name, x, _mesh, *_local_to_global_index_map,
+ _process_variables, _secondary_variables, _process_output);
}
void initialize()
@@ -117,7 +92,8 @@ public:
computeSparsityPattern();
#endif
- createAssemblers(*_local_to_global_index_map, _mesh, _integration_order);
+ initializeConcreteProcess(*_local_to_global_index_map, _mesh,
+ _integration_order);
DBUG("Initialize boundary conditions.");
for (ProcessVariable& pv : _process_variables)
@@ -203,7 +179,7 @@ public:
private:
/// Process specific initialization called by initialize().
- virtual void createAssemblers(
+ virtual void initializeConcreteProcess(
AssemblerLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh,
unsigned const integration_order) = 0;
@@ -326,7 +302,6 @@ private:
unsigned const _integration_order = 2;
MeshLib::Mesh& _mesh;
- std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
_local_to_global_index_map;
@@ -341,6 +316,12 @@ private:
/// Variables used by this process.
std::vector<std::reference_wrapper<ProcessVariable>> _process_variables;
+
+ ProcessOutput<GlobalVector> _process_output;
+
+protected:
+ std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
+ SecondaryVariableCollection<GlobalVector> _secondary_variables;
};
/// Find process variables in \c variables whose names match the settings under
diff --git a/ProcessLib/ProcessOutput.h b/ProcessLib/ProcessOutput.h
new file mode 100644
index 0000000000000000000000000000000000000000..0e29556db7bf46f6eb3a1c63e900f9ea6032e07c
--- /dev/null
+++ b/ProcessLib/ProcessOutput.h
@@ -0,0 +1,258 @@
+/**
+ * \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_PROCESSOUTPUT_H
+#define PROCESSLIB_PROCESSOUTPUT_H
+
+#include "MeshLib/IO/VtkIO/VtuInterface.h"
+#include "ProcessVariable.h"
+#include "SecondaryVariable.h"
+
+namespace ProcessLib
+{
+
+//! Holds information about which variables to write to output files.
+template <typename GlobalVector>
+struct ProcessOutput final
+{
+ //! Constructs a new instance.
+ ProcessOutput(BaseLib::ConfigTree const& output_config,
+ std::vector<std::reference_wrapper<ProcessVariable>> const&
+ process_variables,
+ SecondaryVariableCollection<GlobalVector> const& secondary_variables)
+ {
+ auto const out_vars = output_config.getConfSubtree("variables");
+
+ for (auto out_var : out_vars.getConfParamList<std::string>("variable"))
+ {
+ if (output_variables.find(out_var) != output_variables.cend())
+ {
+ ERR("output variable `%s' specified more than once.", out_var.c_str());
+ std::abort();
+ }
+
+ auto pred = [&out_var](ProcessVariable const& pv) {
+ return pv.getName() == out_var;
+ };
+
+ // check if out_var is a process variable
+ auto const& pcs_var = std::find_if(
+ process_variables.cbegin(), process_variables.cend(), pred);
+
+ if (pcs_var == process_variables.cend()
+ && !secondary_variables.variableExists(out_var))
+ {
+ ERR("Output variable `%s' is neither a process variable nor a"
+ " secondary variable", out_var.c_str());
+ std::abort();
+ }
+
+ DBUG("adding output variable `%s'", out_var.c_str());
+ output_variables.insert(out_var);
+ }
+
+ if (auto out_resid = output_config.getConfParamOptional<bool>(
+ "output_extrapolation_residuals")) {
+ output_residuals = *out_resid;
+ }
+
+ // debug output
+ if (auto const param =
+ output_config.getConfParamOptional<bool>("output_iteration_results"))
+ {
+ DBUG("output_iteration_results: %s", (*param) ? "true" : "false");
+
+ output_iteration_results = *param;
+ }
+ }
+
+ //! All variables that shall be output.
+ std::set<std::string> output_variables;
+
+ //! Tells if also to output extrapolation residuals.
+ bool output_residuals = false;
+
+ bool output_iteration_results = false;
+};
+
+
+//! Writes output to the given \c file_name using the VTU file format.
+template <typename GlobalVector>
+void doProcessOutput(
+ std::string const& file_name,
+ GlobalVector const& x,
+ MeshLib::Mesh& mesh,
+ AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+ std::vector<std::reference_wrapper<ProcessVariable>> const&
+ process_variables,
+ SecondaryVariableCollection<GlobalVector> secondary_variables,
+ ProcessOutput<GlobalVector> const& process_output)
+{
+ DBUG("Process output.");
+
+ // Copy result
+#ifdef USE_PETSC
+ // TODO It is also possible directly to copy the data for single process
+ // variable to a mesh property. It needs a vector of global indices and
+ // some PETSc magic to do so.
+ std::vector<double> x_copy(x.getLocalSize() + x.getGhostSize());
+#else
+ std::vector<double> x_copy(x.size());
+#endif
+ x.copyValues(x_copy);
+
+ auto const& output_variables = process_output.output_variables;
+
+ std::size_t const n_nodes = mesh.getNNodes();
+ int global_component_offset = 0;
+ int global_component_offset_next = 0;
+
+ // primary variables
+ for (ProcessVariable& pv : process_variables)
+ {
+ int const n_components = pv.getNumberOfComponents();
+ global_component_offset = global_component_offset_next;
+ global_component_offset_next += n_components;
+
+ if (output_variables.find(pv.getName()) == output_variables.cend())
+ continue;
+
+ DBUG(" process variable %s", pv.getName().c_str());
+
+ auto& output_data = pv.getOrCreateMeshProperty();
+
+ for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
+ {
+ MeshLib::Location const l(mesh.getID(),
+ MeshLib::MeshItemType::Node, node_id);
+ // TODO extend component ids to multiple process variables.
+ for (int component_id = 0; component_id < n_components;
+ ++component_id)
+ {
+ auto const global_component_id = global_component_offset + component_id;
+ auto const index =
+ dof_table.getLocalIndex(
+ l, global_component_id, x.getRangeBegin(),
+ x.getRangeEnd());
+
+ output_data[node_id * n_components + component_id] =
+ x_copy[index];
+ }
+ }
+ }
+
+#ifndef USE_PETSC
+ // the following section is for the output of secondary variables
+
+ auto count_mesh_items = [](
+ MeshLib::Mesh const& mesh, MeshLib::MeshItemType type) -> std::size_t
+ {
+ switch (type) {
+ case MeshLib::MeshItemType::Cell: return mesh.getNElements();
+ case MeshLib::MeshItemType::Node: return mesh.getNNodes();
+ default: break; // avoid compiler warning
+ }
+ return 0;
+ };
+
+ auto get_or_create_mesh_property = [&mesh, &count_mesh_items](
+ std::string const& property_name, MeshLib::MeshItemType type)
+ {
+ // Get or create a property vector for results.
+ boost::optional<MeshLib::PropertyVector<double>&> result;
+
+ auto const N = count_mesh_items(mesh, type);
+
+ if (mesh.getProperties().hasPropertyVector(property_name))
+ {
+ result = mesh.getProperties().template
+ getPropertyVector<double>(property_name);
+ }
+ else
+ {
+ result = mesh.getProperties().template
+ createNewPropertyVector<double>(property_name, type);
+ result->resize(N);
+ }
+ assert(result && result->size() == N);
+
+ return result;
+ };
+
+ auto add_secondary_var = [&](SecondaryVariable<GlobalVector> const& var,
+ std::string const& output_name)
+ {
+ assert(var.n_components == 1); // TODO implement other cases
+
+ {
+ DBUG(" secondary variable %s", output_name.c_str());
+
+ auto result = get_or_create_mesh_property(
+ output_name, MeshLib::MeshItemType::Node);
+ assert(result->size() == mesh.getNNodes());
+
+ std::unique_ptr<GlobalVector> result_cache;
+ auto const& nodal_values =
+ var.fcts.eval_field(x, dof_table, result_cache);
+
+ // Copy result
+ for (std::size_t i = 0; i < mesh.getNNodes(); ++i)
+ {
+ assert(!std::isnan(nodal_values[i]));
+ (*result)[i] = nodal_values[i];
+ }
+ }
+
+ if (process_output.output_residuals && var.fcts.eval_residuals)
+ {
+ DBUG(" secondary variable %s residual", output_name.c_str());
+ auto const& property_name_res = output_name + "_residual";
+
+ auto result = get_or_create_mesh_property(
+ property_name_res, MeshLib::MeshItemType::Cell);
+ assert(result->size() == mesh.getNElements());
+
+ std::unique_ptr<GlobalVector> result_cache;
+ auto const& residuals =
+ var.fcts.eval_residuals(x, dof_table, result_cache);
+
+ // Copy result
+ for (std::size_t i = 0; i < mesh.getNElements(); ++i)
+ {
+ assert(!std::isnan(residuals[i]));
+ (*result)[i] = residuals[i];
+ }
+ }
+ };
+
+ for (auto const& var : secondary_variables)
+ {
+ auto const var_name = secondary_variables.getMappedName(var.name);
+ if (output_variables.find(var_name)
+ != output_variables.cend())
+ {
+ add_secondary_var(var, var_name);
+ }
+ }
+
+ // secondary variables output end
+#else
+ (void) secondary_variables;
+#endif // USE_PETSC
+
+ // Write output file
+ DBUG("Writing output to \'%s\'.", file_name.c_str());
+ MeshLib::IO::VtuInterface vtu_interface(&mesh, vtkXMLWriter::Binary, true);
+ vtu_interface.writeToFile(file_name);
+}
+
+} // ProcessLib
+
+
+#endif // PROCESSLIB_PROCESSOUTPUT_H
diff --git a/ProcessLib/SecondaryVariable.h b/ProcessLib/SecondaryVariable.h
new file mode 100644
index 0000000000000000000000000000000000000000..7e3cbfbef0c2ccd7968baa9a8444e5a44f53a172
--- /dev/null
+++ b/ProcessLib/SecondaryVariable.h
@@ -0,0 +1,237 @@
+/**
+ * \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_SECONDARY_VARIABLE_H
+#define PROCESSLIB_SECONDARY_VARIABLE_H
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/uniqueInsert.h"
+#include "NumLib/Extrapolation/Extrapolator.h"
+
+namespace AssemblerLib { class LocalToGlobalIndexMap; }
+
+namespace ProcessLib
+{
+
+//! Holder for function objects that compute secondary variables,
+//! and (optionally) also the residuals (e.g., in case of extrapolation)
+template<typename GlobalVector>
+struct SecondaryVariableFunctions final
+{
+ /*! Type of functions used.
+ *
+ * \note The argument \c dof_table is the d.o.f. table of the process, i.e.
+ * it possibly contains information about several process variables.
+ *
+ * \remark The \c result_cache can be used to store the \c GlobalVector if it
+ * is computed on-the-fly. Then a reference to the result cache can be returned.
+ * Otherwise the \c Function must return a reference to a \c GlobalVector that
+ * is stored somewhere else.
+ */
+ using Function = std::function<GlobalVector const&(
+ GlobalVector const& x,
+ AssemblerLib::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_))
+ {
+ // Used to detect nasty implicit conversions.
+ static_assert(std::is_same<GlobalVector const&,
+ typename std::result_of<F1(
+ GlobalVector const&, AssemblerLib::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&, AssemblerLib::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_))
+ {
+ // Used to detect nasty implicit conversions.
+ static_assert(std::is_same<GlobalVector const&,
+ typename std::result_of<F1(
+ GlobalVector const&, AssemblerLib::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;
+};
+
+//! Stores information about a specific secondary variable
+template<typename GlobalVector>
+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<GlobalVector> fcts;
+};
+
+//! Handles configuration of several secondary variables from the project file.
+template<typename GlobalVector>
+class SecondaryVariableCollection final
+{
+public:
+ /*! Constructs new instance.
+ *
+ * \param config Configuration settings.
+ * \param tag_names Possible tag names that contain information about specific
+ * secondary variables.
+ *
+ * In this method only the mapping from tag names to variables is set up.
+ * Any further information has to be passed via addSecondaryVariable().
+ */
+ SecondaryVariableCollection(
+ boost::optional<BaseLib::ConfigTree> const& config,
+ std::initializer_list<std::string> tag_names)
+ {
+ if (!config) return;
+
+ // read which variables are defined in the config
+ for (auto const& tag_name : tag_names) {
+ if (auto var_name = config->getConfParamOptional<std::string>(tag_name))
+ {
+ // TODO check primary vars, too
+ BaseLib::insertIfKeyValueUniqueElseError(
+ _map_tagname_to_varname, tag_name, *var_name,
+ "Secondary variable names must be unique.");
+ }
+ }
+ }
+
+ /*! Tells if a secondary variable with the specified name has been set up.
+ *
+ * \note \c variable_name is not the tag name in the project file!
+ */
+ bool variableExists(std::string const& variable_name) const
+ {
+ auto pred = [&variable_name](std::pair<std::string, std::string> const& p) {
+ return p.second == variable_name;
+ };
+
+ // check if out_var is a secondary variable
+ auto const& var = std::find_if(
+ _map_tagname_to_varname.cbegin(), _map_tagname_to_varname.cend(), pred);
+
+ return var != _map_tagname_to_varname.cend();
+ }
+
+ /*! Set up a secondary variable.
+ *
+ * \param tag_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
+ * Only variables requested by the user in the project file will be configured.
+ * All other variables are silently ignored.
+ */
+ void addSecondaryVariable(
+ std::string const& tag_name, const unsigned num_components,
+ SecondaryVariableFunctions<GlobalVector>&& fcts)
+ {
+ auto it = _map_tagname_to_varname.find(tag_name);
+
+ // get user-supplied var_name for the given tag_name
+ if (it != _map_tagname_to_varname.end()) {
+ auto const& var_name = it->first;
+ // TODO make sure the same variable is not pushed twice
+ _secondary_variables.push_back(
+ {var_name, num_components, std::move(fcts)});
+ }
+ }
+
+ std::string const& getMappedName(std::string const& tag_name)
+ {
+ return _map_tagname_to_varname[tag_name];
+ }
+
+ //! Returns an iterator to the first secondary variable.
+ typename std::vector<SecondaryVariable<GlobalVector>>::const_iterator
+ begin() const
+ {
+ return _secondary_variables.begin();
+ }
+
+ //! Returns an iterator past the last secondary variable.
+ typename std::vector<SecondaryVariable<GlobalVector>>::const_iterator
+ end() const
+ {
+ return _secondary_variables.end();
+ }
+
+private:
+ //! Maps project file tag names to secondary variable names.
+ std::map<std::string, std::string> _map_tagname_to_varname;
+
+ //! Collection of all configured secondary variables.
+ std::vector<SecondaryVariable<GlobalVector>> _secondary_variables;
+};
+
+
+//! Creates an object that computes a secondary variable via extrapolation
+//! of integration point values.
+template<typename GlobalVector, typename PropertyEnum, typename LocalAssembler>
+SecondaryVariableFunctions<GlobalVector>
+makeExtrapolator(PropertyEnum const property,
+ NumLib::Extrapolator<GlobalVector, PropertyEnum, LocalAssembler>&
+ extrapolator,
+ typename NumLib::Extrapolator<GlobalVector, PropertyEnum,
+ LocalAssembler>::LocalAssemblers const& local_assemblers)
+{
+ static_assert(std::is_base_of<
+ NumLib::Extrapolatable<GlobalVector, PropertyEnum>, LocalAssembler>::value,
+ "The passed local assembler type (i.e. the local assembler interface) must"
+ " derive from NumLib::Extrapolatable<>.");
+
+ auto const eval_field = [property, &extrapolator, &local_assemblers](
+ GlobalVector const& /*x*/,
+ AssemblerLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::unique_ptr<GlobalVector>& /*result_cache*/
+ ) -> GlobalVector const&
+ {
+ extrapolator.extrapolate(local_assemblers, property);
+ return extrapolator.getNodalValues();
+ };
+
+ auto const eval_residuals = [property, &extrapolator, &local_assemblers](
+ GlobalVector const& /*x*/,
+ AssemblerLib::LocalToGlobalIndexMap const& /*dof_table*/,
+ std::unique_ptr<GlobalVector>& /*result_cache*/
+ ) -> GlobalVector const&
+ {
+ extrapolator.calculateResiduals(local_assemblers, property);
+ return extrapolator.getElementResiduals();
+ };
+ return { eval_field, eval_residuals };
+}
+
+} // namespace ProcessLib
+
+#endif // PROCESSLIB_SECONDARY_VARIABLE_H
diff --git a/Tests/Data b/Tests/Data
index b45770086122b3264fc47584e3cba7d56698ac54..10fa58f86f4c7c87914e6f94c42751ed8c184392 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit b45770086122b3264fc47584e3cba7d56698ac54
+Subproject commit 10fa58f86f4c7c87914e6f94c42751ed8c184392