From 24cdce4b0a2b9219aa6a72b2503dbb63ec31106d Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Wed, 11 May 2016 17:03:08 -0400
Subject: [PATCH] [PL] added class for output and secondary vars
---
ProcessLib/ProcessOutput.h | 258 +++++++++++++++++++++++++++++++++
ProcessLib/SecondaryVariable.h | 237 ++++++++++++++++++++++++++++++
2 files changed, 495 insertions(+)
create mode 100644 ProcessLib/ProcessOutput.h
create mode 100644 ProcessLib/SecondaryVariable.h
diff --git a/ProcessLib/ProcessOutput.h b/ProcessLib/ProcessOutput.h
new file mode 100644
index 00000000000..0e29556db7b
--- /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 00000000000..7e3cbfbef0c
--- /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
--
GitLab