From 4facfada7039f59b550a30ed632737135798631e Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 16 Aug 2016 16:48:00 +0200
Subject: [PATCH] [PL] IC essentially is a Parameter
---
ProcessLib/InitialCondition.cpp | 103 --------------------------------
ProcessLib/InitialCondition.h | 77 +++---------------------
ProcessLib/Process.cpp | 66 ++++++++++----------
ProcessLib/Process.h | 7 ---
ProcessLib/ProcessVariable.cpp | 47 +++++----------
ProcessLib/ProcessVariable.h | 31 +++-------
ProcessLib/Utils/ProcessUtils.h | 70 +++++++++++++---------
7 files changed, 106 insertions(+), 295 deletions(-)
delete mode 100644 ProcessLib/InitialCondition.cpp
diff --git a/ProcessLib/InitialCondition.cpp b/ProcessLib/InitialCondition.cpp
deleted file mode 100644
index 9eb5bcf5e1e..00000000000
--- a/ProcessLib/InitialCondition.cpp
+++ /dev/null
@@ -1,103 +0,0 @@
-/**
- * \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 "InitialCondition.h"
-
-#include <boost/optional.hpp>
-#include <logog/include/logog.hpp>
-
-#include "BaseLib/ConfigTree.h"
-#include "BaseLib/Error.h"
-
-#include "MathLib/Point3d.h"
-#include "MeshLib/Elements/Element.h"
-#include "MeshLib/Mesh.h"
-
-
-namespace ProcessLib
-{
-std::unique_ptr<InitialCondition> createUniformInitialCondition(
- BaseLib::ConfigTree const& config, int const n_components)
-{
- //! \ogs_file_param{initial_condition__type}
- config.checkConfigParameter("type", "Uniform");
-
- // Optional case for single-component variables where the value can be used.
- // If the 'value' tag is found, use it and return. Otherwise continue with
- // then required tag 'values'.
- if (n_components == 1)
- {
- //! \ogs_file_param{initial_condition__Uniform__value}
- auto const value = config.getConfigParameterOptional<double>("value");
- if (value)
- {
- DBUG("Using value %g for the initial condition.", *value);
- return std::unique_ptr<InitialCondition>{
- new UniformInitialCondition{std::vector<double>{*value}}};
- }
- }
-
- std::vector<double> const values =
- //! \ogs_file_param{initial_condition__Uniform__values}
- config.getConfigParameter<std::vector<double>>("values");
- if (values.size() != static_cast<std::size_t>(n_components))
- {
- OGS_FATAL(
- "Number of values for the initial condition is different from the "
- "number of components.");
- }
-
- DBUG("Using following values for the initial condition:");
- for (double const v : values)
- {
- (void)v; // unused value if building w/o DBUG output.
- DBUG("\t%g", v);
- }
-
- return std::unique_ptr<InitialCondition>{
- new UniformInitialCondition{values}};
-}
-
-std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition(
- BaseLib::ConfigTree const& config,
- MeshLib::Mesh const& mesh,
- int const n_components)
-{
- //! \ogs_file_param{initial_condition__type}
- config.checkConfigParameter("type", "MeshProperty");
-
- //! \ogs_file_param{initial_condition__MeshProperty__field_name}
- auto field_name = config.getConfigParameter<std::string>("field_name");
- DBUG("Using field_name %s", field_name.c_str());
-
- if (!mesh.getProperties().hasPropertyVector(field_name))
- {
- OGS_FATAL("The required property %s does not exists in the mesh.",
- field_name.c_str());
- }
- auto const& property =
- mesh.getProperties().template getPropertyVector<double>(field_name);
- if (!property)
- {
- OGS_FATAL("The required property %s is not of the requested type.",
- field_name.c_str());
- }
-
- if (property->getNumberOfComponents() !=
- static_cast<std::size_t>(n_components))
- {
- OGS_FATAL("The required property %s has different number of components %d, "
- "expected %d.",
- field_name.c_str(), property->getNumberOfComponents(), n_components);
- }
- return std::unique_ptr<InitialCondition>(
- new MeshPropertyInitialCondition(*property));
-}
-
-} // namespace ProcessLib
diff --git a/ProcessLib/InitialCondition.h b/ProcessLib/InitialCondition.h
index a8e8fc1d2b9..60a89bfcca9 100644
--- a/ProcessLib/InitialCondition.h
+++ b/ProcessLib/InitialCondition.h
@@ -10,89 +10,26 @@
#ifndef PROCESS_LIB_INITIAL_CONDITION_H_
#define PROCESS_LIB_INITIAL_CONDITION_H_
-#include <cassert>
-#include "BaseLib/ConfigTree.h"
-#include "MeshLib/PropertyVector.h"
-
-namespace BaseLib
-{
-class ConfigTree;
-}
-
-namespace MeshLib
-{
-template <typename>
-class PropertyVector;
-class Mesh;
-}
+#include "Parameter.h"
namespace ProcessLib
{
-/// The InitialCondition is a base class for spatial distributions of values
-/// defined on mesh nodes.
+// TODO document
class InitialCondition
{
public:
- virtual ~InitialCondition() = default;
- virtual double getValue(std::size_t const node_id,
- int const component_id) const = 0;
-};
-
-/// Uniform value initial condition
-class UniformInitialCondition : public InitialCondition
-{
-public:
- explicit UniformInitialCondition(std::vector<double> const& values)
- : _values(values)
- {
- }
- /// Returns a value for given node and component.
- virtual double getValue(std::size_t const /*node_id*/,
- int const component_id) const override
- {
- return _values[component_id];
- }
-
-private:
- std::vector<double> const _values;
-};
-
-/// Construct a UniformInitialCondition from configuration.
-/// The initial condition will expect a correct number of components in the
-/// configuration, which should be the same as in the corresponding process
-/// variable.
-std::unique_ptr<InitialCondition> createUniformInitialCondition(
- BaseLib::ConfigTree const& config, int const n_components);
-
-/// Distribution of values given by a mesh property defined on nodes.
-class MeshPropertyInitialCondition : public InitialCondition
-{
-public:
- MeshPropertyInitialCondition(
- MeshLib::PropertyVector<double> const& property)
- : _property(property)
- {
- assert(_property.getMeshItemType() == MeshLib::MeshItemType::Node);
- }
+ InitialCondition(Parameter<double> const& param) : _param(param) {}
- virtual double getValue(std::size_t const node_id,
- int const component_id) const override
+ std::vector<double> const& getTuple(
+ double const t, SpatialPosition const& pos) const
{
- return _property.getComponent(node_id, component_id);
+ return _param.getTuple(t, pos);
}
private:
- MeshLib::PropertyVector<double> const& _property;
+ Parameter<double> const& _param;
};
-/// Construct a MeshPropertyInitialCondition from configuration.
-/// The initial condition will expect a correct number of components in the
-/// configuration, which should be the same as in the corresponding process
-/// variable.
-std::unique_ptr<InitialCondition> createMeshPropertyInitialCondition(
- BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh,
- int const n_components);
-
} // namespace ProcessLib
#endif // PROCESS_LIB_INITIAL_CONDITION_H_
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 4a82d136519..84d1c17eb88 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -72,15 +72,47 @@ void Process::initialize()
void Process::setInitialConditions(GlobalVector& x)
{
DBUG("Set initial conditions.");
+ std::size_t const n_nodes = _mesh.getNumberOfNodes();
+
+ SpatialPosition pos;
+
for (int variable_id = 0;
variable_id < static_cast<int>(_process_variables.size());
++variable_id)
{
ProcessVariable& pv = _process_variables[variable_id];
- for (int component_id = 0; component_id < pv.getNumberOfComponents();
- ++component_id)
+ auto const* ic = pv.getInitialCondition();
+ if (!ic)
+ continue;
+
+ auto const num_comp = pv.getNumberOfComponents();
+
+ for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
{
- setInitialConditions(pv, variable_id, component_id, x);
+ MeshLib::Location const l(_mesh.getID(),
+ MeshLib::MeshItemType::Node, node_id);
+
+ pos.setNodeID(node_id);
+ auto const& tup = ic->getTuple(0.0, pos); // 0.0 is t!
+
+ for (int comp_id = 0; comp_id < num_comp; ++comp_id)
+ {
+ auto global_index =
+ std::abs(_local_to_global_index_map->getGlobalIndex(
+ l, variable_id, comp_id));
+#ifdef USE_PETSC
+ // The global indices of the ghost entries of the global matrix
+ // or the global vectors need to be set as negative values for
+ // equation assembly, however the global indices start from
+ // zero. Therefore, any ghost entry with zero index is assigned
+ // an negative value of the vector size or the matrix dimension.
+ // To assign the initial value for the ghost entries, the
+ // negative indices of the ghost entries are restored to zero.
+ if (global_index == x.size())
+ global_index = 0;
+#endif
+ x.set(global_index, tup[comp_id]);
+ }
}
}
}
@@ -224,34 +256,6 @@ void Process::finishNamedFunctionsInitialization()
}
}
-void Process::setInitialConditions(ProcessVariable const& variable,
- int const variable_id,
- int const component_id,
- GlobalVector& x)
-{
- std::size_t const n_nodes = _mesh.getNumberOfNodes();
- for (std::size_t node_id = 0; node_id < n_nodes; ++node_id)
- {
- MeshLib::Location const l(_mesh.getID(), MeshLib::MeshItemType::Node,
- node_id);
- auto global_index = std::abs(_local_to_global_index_map->getGlobalIndex(
- l, variable_id, component_id));
-#ifdef USE_PETSC
- // The global indices of the ghost entries of the global matrix or
- // the global vectors need to be set as negative values for equation
- // assembly, however the global indices start from zero. Therefore,
- // any ghost entry with zero index is assigned an negative value of
- // the vector size or the matrix dimension. To assign the initial
- // value for the ghost entries, the negative indices of the ghost
- // entries are restored to zero.
- if (global_index == x.size())
- global_index = 0;
-#endif
- x.set(global_index,
- variable.getInitialConditionValue(node_id, component_id));
- }
-}
-
void Process::computeSparsityPattern()
{
_sparsity_pattern =
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 2d9595abfdf..7879580d639 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -148,13 +148,6 @@ private:
/// each of the named functions.
void finishNamedFunctionsInitialization();
- /// Sets the initial condition values in the solution vector x for a given
- /// process variable and component.
- void setInitialConditions(ProcessVariable const& variable,
- int const variable_id,
- int const component_id,
- GlobalVector& x);
-
/// Computes and stores global matrix' sparsity pattern from given
/// DOF-table.
void computeSparsityPattern();
diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp
index 618c2857534..40d2645adcf 100644
--- a/ProcessLib/ProcessVariable.cpp
+++ b/ProcessLib/ProcessVariable.cpp
@@ -10,48 +10,33 @@
#include "ProcessVariable.h"
#include <utility>
-
#include <logog/include/logog.hpp>
#include "GeoLib/GEOObjects.h"
#include "MeshLib/Mesh.h"
-
-#include "BaseLib/ConfigTree.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
namespace ProcessLib
{
-ProcessVariable::ProcessVariable(BaseLib::ConfigTree const& config,
- MeshLib::Mesh& mesh,
- GeoLib::GEOObjects const& geometries)
- :
- //! \ogs_file_param{prj__process_variables__process_variable__name}
- _name(config.getConfigParameter<std::string>("name")),
- _mesh(mesh),
- //! \ogs_file_param{prj__process_variables__process_variable__components}
- _n_components(config.getConfigParameter<int>("components"))
+ProcessVariable::ProcessVariable(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
+ GeoLib::GEOObjects const& geometries,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters)
+ : //! \ogs_file_param{prj__process_variables__process_variable__name}
+ _name(config.getConfigParameter<std::string>("name")),
+ _mesh(mesh),
+ //! \ogs_file_param{prj__process_variables__process_variable__components}
+ _n_components(config.getConfigParameter<int>("components"))
{
- DBUG("Constructing process variable %s", this->_name.c_str());
+ DBUG("Constructing process variable %s", _name.c_str());
// Initial condition
- //! \ogs_file_param{prj__process_variables__process_variable__initial_condition}
- if (auto ic_config = config.getConfigSubtreeOptional("initial_condition"))
+ if (auto ic_name =
+ //! \ogs_file_param{prj__process_variables__process_variable__initial_condition}
+ config.getConfigParameterOptional<std::string>("initial_condition"))
{
- //! \ogs_file_param{initial_condition__type}
- auto const type = ic_config->peekConfigParameter<std::string>("type");
- if (type == "Uniform")
- {
- _initial_condition =
- createUniformInitialCondition(*ic_config, _n_components);
- }
- else if (type == "MeshProperty")
- {
- _initial_condition =
- createMeshPropertyInitialCondition(*ic_config, _mesh, _n_components);
- }
- else
- {
- ERR("Unknown type of the initial condition.");
- }
+ _initial_condition.reset(new InitialCondition(
+ findParameter<double>(*ic_name, parameters, _n_components)));
}
else
{
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index 463c9d1a271..ca2828e2acd 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -10,32 +10,14 @@
#ifndef PROCESS_LIB_PROCESS_VARIABLE_H_
#define PROCESS_LIB_PROCESS_VARIABLE_H_
-
#include "InitialCondition.h"
#include "ProcessLib/BoundaryCondition/BoundaryCondition.h"
#include "ProcessLib/BoundaryCondition/BoundaryConditionConfig.h"
-namespace MeshGeoToolsLib
-{
-class MeshNodeSearcher;
-class BoundaryElementsSearcher;
-}
-
namespace MeshLib
{
class Mesh;
-}
-
-namespace GeoLib
-{
-class GEOObjects;
-}
-
-namespace ProcessLib
-{
-class NeumannBcConfig;
-class InitialCondition;
-class UniformDirichletBoundaryCondition;
+template <typename T> class PropertyVector;
}
namespace ProcessLib
@@ -45,8 +27,10 @@ namespace ProcessLib
class ProcessVariable
{
public:
- ProcessVariable(BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
- GeoLib::GEOObjects const& geometries);
+ ProcessVariable(
+ BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh,
+ GeoLib::GEOObjects const& geometries,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters);
ProcessVariable(ProcessVariable&&);
@@ -63,10 +47,9 @@ public:
const int variable_id,
unsigned const integration_order);
- double getInitialConditionValue(std::size_t const node_id,
- int const component_id) const
+ InitialCondition const* getInitialCondition() const
{
- return _initial_condition->getValue(node_id, component_id);
+ return _initial_condition.get();
}
// Get or create a property vector for results.
diff --git a/ProcessLib/Utils/ProcessUtils.h b/ProcessLib/Utils/ProcessUtils.h
index de0e73ffacd..57062c652fa 100644
--- a/ProcessLib/Utils/ProcessUtils.h
+++ b/ProcessLib/Utils/ProcessUtils.h
@@ -43,42 +43,26 @@ std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables(
std::initializer_list<std::string>
tag_names);
-/// Find a parameter of specific type for a name given in the process
-/// configuration under the tag.
-/// The parameter must have the specified number of components.
-/// In the process config a parameter is referenced by a name. For example it
-/// will be looking for a parameter named "K" in the list of parameters
-/// when the tag is "hydraulic_conductivity":
-/// \code
-/// <process>
-/// ...
-/// <hydraulic_conductivity>K</hydraulic_conductivity>
-/// </process>
-/// \endcode
-/// and return a reference to that parameter. Additionally it checks for the
-/// type of the found parameter.
+/// Find a parameter of specific type for a given name.
+///
+/// \see The documentation of the other findParameter() function.
template <typename ParameterDataType>
Parameter<ParameterDataType>& findParameter(
- BaseLib::ConfigTree const& process_config, std::string const& tag,
+ std::string const& parameter_name,
std::vector<std::unique_ptr<ParameterBase>> const& parameters,
unsigned const num_components)
{
- // Find parameter name in process config.
- //! \ogs_file_special
- auto const name = process_config.getConfigParameter<std::string>(tag);
-
// Find corresponding parameter by name.
- auto const parameter_it =
- std::find_if(parameters.cbegin(), parameters.cend(),
- [&name](std::unique_ptr<ParameterBase> const& p) {
- return p->name == name;
- });
+ auto const parameter_it = std::find_if(
+ parameters.cbegin(), parameters.cend(),
+ [¶meter_name](std::unique_ptr<ParameterBase> const& p) {
+ return p->name == parameter_name;
+ });
if (parameter_it == parameters.end()) {
OGS_FATAL(
- "Could not find parameter `%s' in the provided parameters list for "
- "config tag <%s>.",
- name.c_str(), tag.c_str());
+ "Could not find parameter `%s' in the provided parameters list.",
+ parameter_name.c_str());
}
DBUG("Found parameter `%s'.", (*parameter_it)->name.c_str());
@@ -87,18 +71,46 @@ Parameter<ParameterDataType>& findParameter(
dynamic_cast<Parameter<ParameterDataType>*>(parameter_it->get());
if (!parameter) {
OGS_FATAL("The read parameter `%s' is of incompatible type.",
- name.c_str());
+ parameter_name.c_str());
}
if (parameter->getNumberOfComponents() != num_components) {
OGS_FATAL(
"The read parameter `%s' has the wrong number of components (%lu "
"instead of %u).",
- name.c_str(), parameter->getNumberOfComponents(), num_components);
+ parameter_name.c_str(), parameter->getNumberOfComponents(),
+ num_components);
}
return *parameter;
}
+
+/// Find a parameter of specific type for a name given in the process
+/// configuration under the tag.
+/// The parameter must have the specified number of components.
+/// In the process config a parameter is referenced by a name. For example it
+/// will be looking for a parameter named "K" in the list of parameters
+/// when the tag is "hydraulic_conductivity":
+/// \code
+/// <process>
+/// ...
+/// <hydraulic_conductivity>K</hydraulic_conductivity>
+/// </process>
+/// \endcode
+/// and return a reference to that parameter. Additionally it checks for the
+/// type of the found parameter.
+template <typename ParameterDataType>
+Parameter<ParameterDataType>& findParameter(
+ BaseLib::ConfigTree const& process_config, std::string const& tag,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const num_components)
+{
+ // Find parameter name in process config.
+ //! \ogs_file_special
+ auto const name = process_config.getConfigParameter<std::string>(tag);
+
+ return findParameter<ParameterDataType>(name, parameters, num_components);
+}
} // ProcessLib
#endif // PROCESSLIB_UTILS_PROCESSUTILS_H
--
GitLab