diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index 29a2afe89b8601bda5ead02c49f60c7aedee595c..07a48a655155ae4316995fe7a4e82cdb42423412 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -76,12 +76,12 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config, //! \ogs_file_param{prj__curves} parseCurves(project_config.getConfigSubtreeOptional("curves")); - //! \ogs_file_param{prj__process_variables} - parseProcessVariables(project_config.getConfigSubtree("process_variables")); - //! \ogs_file_param{prj__parameters} parseParameters(project_config.getConfigSubtree("parameters")); + //! \ogs_file_param{prj__process_variables} + parseProcessVariables(project_config.getConfigSubtree("process_variables")); + //! \ogs_file_param{prj__processes} parseProcesses(project_config.getConfigSubtree("processes")); @@ -273,7 +273,7 @@ void ProjectData::parseProcessVariables( { // TODO Extend to referenced meshes. _process_variables.emplace_back(var_config, *_mesh_vec[0], - *_geoObjects); + *_geoObjects, _parameters); } } @@ -286,30 +286,8 @@ void ProjectData::parseParameters(BaseLib::ConfigTree const& parameters_config) //! \ogs_file_param{prj__parameters__parameter} parameters_config.getConfigSubtreeList("parameter")) { - //! \ogs_file_param{parameter__name} - auto name = parameter_config.getConfigParameter<std::string>("name"); - //! \ogs_file_param{parameter__type} - auto type = parameter_config.peekConfigParameter<std::string>("type"); - - // Create parameter based on the provided type. - if (type == "Constant") - { - INFO("ConstantParameter: %s.", name.c_str()); - _parameters.push_back(createConstParameter(parameter_config)); - _parameters.back()->name = name; - } - else if (type == "MeshProperty") - { - INFO("MeshPropertyParameter: %s", name.c_str()); - _parameters.push_back( - createMeshPropertyParameter(parameter_config, *_mesh_vec[0])); - _parameters.back()->name = name; - } - else - { - OGS_FATAL("Cannot construct property of given type \'%s\'.", - type.c_str()); - } + _parameters.push_back( + ProcessLib::createParameter(parameter_config, _mesh_vec)); } } diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h index 08d9857a1dc9ac17b16be3219651ef0b8e2738c3..ed8006bfd1767858f4b71af1b6b1d8f1cceb373b 100644 --- a/Applications/ApplicationsLib/ProjectData.h +++ b/Applications/ApplicationsLib/ProjectData.h @@ -23,7 +23,7 @@ #include "GeoLib/GEOObjects.h" #include "ProcessLib/Output.h" -#include "ProcessLib/Parameter.h" +#include "ProcessLib/Parameter/Parameter.h" #include "ProcessLib/Process.h" #include "ProcessLib/ProcessVariable.h" diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp index 4eae2e6234e7717d5844329c3cede8b64992212a..cdfdecec8e2740576e9c07d3caf3b9a306669d63 100644 --- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp +++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp @@ -89,7 +89,7 @@ void UncoupledProcessesTimeLoop::setInitialConditions( ode_sys.getMatrixSpecifications())); auto& x0 = *_process_solutions[pcs_idx]; - pcs.setInitialConditions(x0); + pcs.setInitialConditions(t0, x0); MathLib::LinAlg::finalizeAssembly(x0); time_disc.setInitialState(t0, x0); // push IC diff --git a/MeshLib/PropertyVector.h b/MeshLib/PropertyVector.h index df9041f205bc60d5b6bcfb69cbdfaab6829fbec2..e0558b44b330a56b6f4bce25e010c496ce07fd94 100644 --- a/MeshLib/PropertyVector.h +++ b/MeshLib/PropertyVector.h @@ -66,6 +66,15 @@ public: return std::vector<PROP_VAL_TYPE>::size() / _n_components; } + //! Returns the value for the given component stored in the given tuple. + PROP_VAL_TYPE const& getComponent(std::size_t tuple_index, + std::size_t component) const + { + assert(component < _n_components); + assert(tuple_index < getNumberOfTuples()); + return this->operator[](tuple_index* getNumberOfTuples() + component); + } + PropertyVectorBase* clone(std::vector<std::size_t> const& exclude_positions) const { PropertyVector<PROP_VAL_TYPE> *t(new PropertyVector<PROP_VAL_TYPE>(_property_name, diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt index d3529016c7398999db3be5ee31132e22109f131e..d3883b0d7784c4fda37015eb1b3dd52308cd35f5 100644 --- a/ProcessLib/CMakeLists.txt +++ b/ProcessLib/CMakeLists.txt @@ -6,6 +6,9 @@ APPEND_SOURCE_FILES(SOURCES) add_subdirectory(BoundaryCondition) APPEND_SOURCE_FILES(SOURCES BoundaryCondition) +add_subdirectory(Parameter) +APPEND_SOURCE_FILES(SOURCES Parameter) + add_subdirectory(GroundwaterFlow) APPEND_SOURCE_FILES(SOURCES GroundwaterFlow) diff --git a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp index 13f816c95d1535b9b631e1a337b2d98dc3017019..bb7919629c5d01e9f2f18985f6f4462fb4fd83ff 100644 --- a/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp +++ b/ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.cpp @@ -10,6 +10,7 @@ #include "CreateGroundwaterFlowProcess.h" #include "ProcessLib/Utils/ParseSecondaryVariables.h" +#include "ProcessLib/Utils/ProcessUtils.h" #include "GroundwaterFlowProcess.h" #include "GroundwaterFlowProcessData.h" @@ -38,12 +39,11 @@ std::unique_ptr<Process> createGroundwaterFlowProcess( "process_variable"}); // Hydraulic conductivity parameter. - auto& hydraulic_conductivity = findParameter<double, - MeshLib::Element const&>( + auto& hydraulic_conductivity = findParameter<double>( config, //! \ogs_file_param_special{process__GROUNDWATER_FLOW__hydraulic_conductivity} "hydraulic_conductivity", - parameters); + parameters, 1); DBUG("Use \'%s\' as hydraulic conductivity parameter.", hydraulic_conductivity.name.c_str()); diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h index 8f1f8b3be89b2748b60273af3ba362b91b0b4d5e..5adda111bec0df75c9873bb467ab315b17611c84 100644 --- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h +++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h @@ -17,7 +17,7 @@ #include "NumLib/Fem/ShapeMatrixPolicy.h" #include "ProcessLib/LocalAssemblerInterface.h" #include "ProcessLib/LocalAssemblerTraits.h" -#include "ProcessLib/Parameter.h" +#include "ProcessLib/Parameter/Parameter.h" #include "ProcessLib/Utils/InitShapeMatrices.h" #include "GroundwaterFlowProcessData.h" @@ -81,7 +81,7 @@ public: } void assembleConcrete( - double const /*t*/, std::vector<double> const& local_x, + double const t, std::vector<double> const& local_x, NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices, GlobalMatrix& /*M*/, GlobalMatrix& K, GlobalVector& b) override { @@ -91,11 +91,15 @@ public: IntegrationMethod integration_method(_integration_order); unsigned const n_integration_points = integration_method.getNumberOfPoints(); - for (std::size_t ip(0); ip < n_integration_points; ip++) + 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& wp = integration_method.getWeightedPoint(ip); - auto const k = _process_data.hydraulic_conductivity(_element); + auto const k = _process_data.hydraulic_conductivity.getTuple(t, pos).front(); _localA.noalias() += sm.dNdx.transpose() * k * sm.dNdx * sm.detJ * wp.getWeight(); diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcessData.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcessData.h index d148011b142446a696354961f6ee882072d94659..e289867d5ca9481c483dafd17d4083f792f72bd2 100644 --- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcessData.h +++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcessData.h @@ -19,7 +19,7 @@ namespace MeshLib namespace ProcessLib { -template <typename ReturnType, typename... Args> +template <typename T> struct Parameter; namespace GroundwaterFlow @@ -28,9 +28,7 @@ namespace GroundwaterFlow struct GroundwaterFlowProcessData { GroundwaterFlowProcessData( - ProcessLib::Parameter<double, MeshLib::Element const&> const& - hydraulic_conductivity_ - ) + Parameter<double> const& hydraulic_conductivity_) : hydraulic_conductivity(hydraulic_conductivity_) {} @@ -47,7 +45,7 @@ struct GroundwaterFlowProcessData //! Assignments are not needed. void operator=(GroundwaterFlowProcessData&&) = delete; - Parameter<double, MeshLib::Element const&> const& hydraulic_conductivity; + Parameter<double> const& hydraulic_conductivity; }; } // namespace GroundwaterFlow diff --git a/ProcessLib/InitialCondition.cpp b/ProcessLib/InitialCondition.cpp deleted file mode 100644 index 9eb5bcf5e1e60d7724ae26b07907183f93a93412..0000000000000000000000000000000000000000 --- 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 deleted file mode 100644 index 62a50480ad8c74dd34a4c7816b3a44ae72c628b8..0000000000000000000000000000000000000000 --- a/ProcessLib/InitialCondition.h +++ /dev/null @@ -1,99 +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 - * - */ - -#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; -} - -namespace ProcessLib -{ -/// The InitialCondition is a base class for spatial distributions of values -/// defined on mesh nodes. -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); - } - - virtual double getValue(std::size_t const node_id, - int const component_id) const override - { - return _property[node_id * _property.getNumberOfComponents() + - component_id]; - } - -private: - MeshLib::PropertyVector<double> const& _property; -}; - -/// 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/Parameter.cpp b/ProcessLib/Parameter.cpp deleted file mode 100644 index 0ad1e9e89eecfb9da4419d310a575a8fb2d8a3bb..0000000000000000000000000000000000000000 --- a/ProcessLib/Parameter.cpp +++ /dev/null @@ -1,57 +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 "Parameter.h" - -#include <boost/optional.hpp> -#include <logog/include/logog.hpp> - -#include "BaseLib/Error.h" -#include "MeshLib/Elements/Element.h" - -namespace ProcessLib -{ -std::unique_ptr<ParameterBase> createConstParameter( - BaseLib::ConfigTree const& config) -{ - //! \ogs_file_param{parameter__type} - config.checkConfigParameter("type", "Constant"); - //! \ogs_file_param{parameter__Constant__value} - auto value = config.getConfigParameter<double>("value"); - DBUG("Using value %g", value); - - return std::unique_ptr<ParameterBase>(new ConstParameter<double>(value)); -} - -std::unique_ptr<ParameterBase> createMeshPropertyParameter( - BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh) -{ - //! \ogs_file_param{parameter__type} - config.checkConfigParameter("type", "MeshProperty"); - //! \ogs_file_param{parameter__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()); - } - - return std::unique_ptr<ParameterBase>( - new MeshPropertyParameter<double>(*property)); -} -} // namespace ProcessLib diff --git a/ProcessLib/Parameter.h b/ProcessLib/Parameter.h deleted file mode 100644 index 4c52a459ee6e68e9d781138797694aa14735bddc..0000000000000000000000000000000000000000 --- a/ProcessLib/Parameter.h +++ /dev/null @@ -1,94 +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 - * - */ - -#ifndef PROCESS_LIB_PARAMETER_H_ -#define PROCESS_LIB_PARAMETER_H_ - -#include <memory> - -#include <logog/include/logog.hpp> -#include <boost/optional.hpp> - -#include "BaseLib/ConfigTree.h" -#include "MeshLib/Elements/Element.h" - -namespace MeshLib -{ -template <typename T> -class PropertyVector; -} - -namespace ProcessLib -{ -/// Base class for all parameters, not an interface class. This avoids using of -/// void* when storing parameters and convenient destruction. -/// Its property name helps addressing the right parameter. -struct ParameterBase -{ - virtual ~ParameterBase() = default; - - std::string name; -}; - -/// A parameter is representing a value or function of any type. -/// The ReturnType can represent an underlying type of an aggregate type like -/// tuple or matrix (\see tupleSize()). -/// The total number of stored tuples is provided. -template <typename ReturnType, typename... Args> -struct Parameter : public ParameterBase -{ - virtual ~Parameter() = default; - - virtual ReturnType operator()(Args&&... args) const = 0; -}; - -/// Single, constant value parameter. -template <typename ReturnType> -struct ConstParameter final - : public Parameter<ReturnType, MeshLib::Element const&> -{ - ConstParameter(ReturnType value) : _value(value) - { - } - - ReturnType operator()(MeshLib::Element const&) const override - { - return _value; - } - -private: - ReturnType _value; -}; - -std::unique_ptr<ParameterBase> createConstParameter(BaseLib::ConfigTree const& config); - -/// A parameter represented by a mesh property vector. -template <typename ReturnType> -struct MeshPropertyParameter final - : public Parameter<ReturnType, MeshLib::Element const&> -{ - MeshPropertyParameter(MeshLib::PropertyVector<ReturnType> const& property) - : _property(property) - { - } - - ReturnType operator()(MeshLib::Element const& e) const override - { - return _property[e.getID()]; - } - -private: - MeshLib::PropertyVector<ReturnType> const& _property; -}; - -std::unique_ptr<ParameterBase> createMeshPropertyParameter(BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh); - -} // namespace ProcessLib - -#endif // PROCESS_LIB_PARAMETER_H_ diff --git a/ProcessLib/Parameter/CMakeLists.txt b/ProcessLib/Parameter/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/ProcessLib/Parameter/ConstantParameter.cpp b/ProcessLib/Parameter/ConstantParameter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ab3e051c746920f83b79e440d9b57f8eb35cb8ad --- /dev/null +++ b/ProcessLib/Parameter/ConstantParameter.cpp @@ -0,0 +1,28 @@ +/** + * \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 "ConstantParameter.h" +#include <logog/include/logog.hpp> +#include "BaseLib/ConfigTree.h" + +namespace ProcessLib +{ +std::unique_ptr<ParameterBase> createConstantParameter( + BaseLib::ConfigTree const& config) +{ + //! \ogs_file_param{parameter__type} + config.checkConfigParameter("type", "Constant"); + //! \ogs_file_param{parameter__Constant__value} + auto const value = config.getConfigParameter<double>("value"); + DBUG("Using value %g", value); + + return std::unique_ptr<ParameterBase>(new ConstantParameter<double>(value)); +} + +} // ProcessLib diff --git a/ProcessLib/Parameter/ConstantParameter.h b/ProcessLib/Parameter/ConstantParameter.h new file mode 100644 index 0000000000000000000000000000000000000000..93007466926312509a7835f025ab3189a052bb74 --- /dev/null +++ b/ProcessLib/Parameter/ConstantParameter.h @@ -0,0 +1,40 @@ +/** + * \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_CONSTANTPARAMETER_H +#define PROCESSLIB_CONSTANTPARAMETER_H + +#include "Parameter.h" + +namespace ProcessLib +{ +/// Single, constant value parameter. +template <typename T> +struct ConstantParameter final : public Parameter<T> { + ConstantParameter(T const& value) : _value({value}) {} + + // TODO allow for different sizes + unsigned getNumberOfComponents() const { return 1; } + + std::vector<T> const& getTuple( + double const /*t*/, SpatialPosition const& /*pos*/) const override + { + return _value; + } + +private: + std::vector<T> const _value; +}; + +std::unique_ptr<ParameterBase> createConstantParameter( + BaseLib::ConfigTree const& config); + +} // ProcessLib + +#endif // PROCESSLIB_CONSTANTPARAMETER_H diff --git a/ProcessLib/Parameter/MeshElementParameter.cpp b/ProcessLib/Parameter/MeshElementParameter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d21ffdeb9dfa110ff40d1ef74a6392d6c29f469b --- /dev/null +++ b/ProcessLib/Parameter/MeshElementParameter.cpp @@ -0,0 +1,48 @@ +/** + * \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 "MeshElementParameter.h" +#include "BaseLib/ConfigTree.h" +#include "BaseLib/Error.h" +#include "MeshLib/Mesh.h" + +namespace ProcessLib +{ +std::unique_ptr<ParameterBase> createMeshElementParameter( + BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh) +{ + //! \ogs_file_param{parameter__type} + config.checkConfigParameter("type", "MeshElement"); + //! \ogs_file_param{parameter__MeshElement__field_name} + auto const 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()); + } + + // TODO other data types than only double + auto const& property = + mesh.getProperties().getPropertyVector<double>(field_name); + if (!property) { + OGS_FATAL("The mesh property `%s' is not of the requested type.", + field_name.c_str()); + } + + if (property->getMeshItemType() != MeshLib::MeshItemType::Cell) { + OGS_FATAL("The mesh property `%s' is not an element property.", + field_name.c_str()); + } + + return std::unique_ptr<ParameterBase>( + new MeshElementParameter<double>(*property)); +} + +} // ProcessLib diff --git a/ProcessLib/Parameter/MeshElementParameter.h b/ProcessLib/Parameter/MeshElementParameter.h new file mode 100644 index 0000000000000000000000000000000000000000..298065a64c7a8cf3599b64c7b3d787a821d7007a --- /dev/null +++ b/ProcessLib/Parameter/MeshElementParameter.h @@ -0,0 +1,59 @@ +/** + * \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_MESHELEMENTPARAMETER_H +#define PROCESSLIB_MESHELEMENTPARAMETER_H + +#include "Parameter.h" + +namespace MeshLib +{ +template <typename T> +class PropertyVector; +} // MeshLib + +namespace ProcessLib +{ +/// A parameter represented by a mesh property vector. +template <typename T> +struct MeshElementParameter final : public Parameter<T> { + MeshElementParameter(MeshLib::PropertyVector<T> const& property) + : _property(property) + , _cache(_property.getNumberOfComponents()) + { + } + + unsigned getNumberOfComponents() const override + { + return _property.getNumberOfComponents(); + } + + std::vector<T> const& getTuple(double const /*t*/, + SpatialPosition const& pos) const override + { + auto const e = pos.getElementID(); + assert(e); + auto const num_comp = _property.getNumberOfComponents(); + for (std::size_t c=0; c<num_comp; ++c) { + _cache[c] = _property.getComponent(*e, c); + } + return _cache; + } + +private: + MeshLib::PropertyVector<T> const& _property; + mutable std::vector<double> _cache; +}; + +std::unique_ptr<ParameterBase> createMeshElementParameter( + BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh); + +} // ProcessLib + +#endif // PROCESSLIB_MESHELEMENTPARAMETER_H diff --git a/ProcessLib/Parameter/MeshNodeParameter.cpp b/ProcessLib/Parameter/MeshNodeParameter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..af986ad603c0471b061dd259971c9090e93d5ccc --- /dev/null +++ b/ProcessLib/Parameter/MeshNodeParameter.cpp @@ -0,0 +1,48 @@ +/** + * \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 "MeshNodeParameter.h" +#include "BaseLib/ConfigTree.h" +#include "BaseLib/Error.h" +#include "MeshLib/Mesh.h" + +namespace ProcessLib +{ +std::unique_ptr<ParameterBase> createMeshNodeParameter( + BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh) +{ + //! \ogs_file_param{parameter__type} + config.checkConfigParameter("type", "MeshNode"); + //! \ogs_file_param{parameter__MeshNode__field_name} + auto const 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()); + } + + // TODO other data types than only double + auto const& property = + mesh.getProperties().getPropertyVector<double>(field_name); + if (!property) { + OGS_FATAL("The mesh property `%s' is not of the requested type.", + field_name.c_str()); + } + + if (property->getMeshItemType() != MeshLib::MeshItemType::Node) { + OGS_FATAL("The mesh property `%s' is not a nodal property.", + field_name.c_str()); + } + + return std::unique_ptr<ParameterBase>( + new MeshNodeParameter<double>(*property)); +} + +} // ProcessLib diff --git a/ProcessLib/Parameter/MeshNodeParameter.h b/ProcessLib/Parameter/MeshNodeParameter.h new file mode 100644 index 0000000000000000000000000000000000000000..2ebecc428de4e0c1d82f7c87ee48e2b60714477e --- /dev/null +++ b/ProcessLib/Parameter/MeshNodeParameter.h @@ -0,0 +1,59 @@ +/** + * \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_MESHNODEPARAMETER_H +#define PROCESSLIB_MESHNODEPARAMETER_H + +#include "Parameter.h" + +namespace MeshLib +{ +template <typename T> +class PropertyVector; +} // MeshLib + +namespace ProcessLib +{ +/// A parameter represented by a mesh property vector. +template <typename T> +struct MeshNodeParameter final : public Parameter<T> { + MeshNodeParameter(MeshLib::PropertyVector<T> const& property) + : _property(property) + , _cache(_property.getNumberOfComponents()) + { + } + + unsigned getNumberOfComponents() const override + { + return _property.getNumberOfComponents(); + } + + std::vector<T> const& getTuple(double const /*t*/, + SpatialPosition const& pos) const override + { + auto const n = pos.getNodeID(); + assert(n); + auto const num_comp = _property.getNumberOfComponents(); + for (std::size_t c=0; c<num_comp; ++c) { + _cache[c] = _property.getComponent(*n, c); + } + return _cache; + } + +private: + MeshLib::PropertyVector<T> const& _property; + mutable std::vector<double> _cache; +}; + +std::unique_ptr<ParameterBase> createMeshNodeParameter( + BaseLib::ConfigTree const& config, MeshLib::Mesh const& mesh); + +} // ProcessLib + +#endif // PROCESSLIB_MESHNODEPARAMETER_H diff --git a/ProcessLib/Parameter/Parameter.cpp b/ProcessLib/Parameter/Parameter.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8655e90182882c6bbcd61910102151baea26876d --- /dev/null +++ b/ProcessLib/Parameter/Parameter.cpp @@ -0,0 +1,59 @@ +/** + * \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 "Parameter.h" +#include "BaseLib/ConfigTree.h" +#include "BaseLib/Error.h" + +#include "ConstantParameter.h" +#include "MeshElementParameter.h" +#include "MeshNodeParameter.h" + +namespace ProcessLib +{ + +std::unique_ptr<ParameterBase> createParameter( + BaseLib::ConfigTree const& config, std::vector<MeshLib::Mesh*> const& meshes) +{ + + //! \ogs_file_param{parameter__name} + auto const name = config.getConfigParameter<std::string>("name"); + //! \ogs_file_param{parameter__type} + auto const type = config.peekConfigParameter<std::string>("type"); + + // Create parameter based on the provided type. + if (type == "Constant") + { + INFO("ConstantParameter: %s", name.c_str()); + auto param = createConstantParameter(config); + param->name = name; + return param; + } + else if (type == "MeshElement") + { + INFO("MeshElementParameter: %s", name.c_str()); + auto param = createMeshElementParameter(config, *meshes.front()); + param->name = name; + return param; + } + else if (type == "MeshNode") + { + INFO("MeshElementParameter: %s", name.c_str()); + auto param = createMeshNodeParameter(config, *meshes.front()); + param->name = name; + return param; + } + else + { + OGS_FATAL("Cannot construct a parameter of given type \'%s\'.", + type.c_str()); + } +} + +} // ProcessLib diff --git a/ProcessLib/Parameter/Parameter.h b/ProcessLib/Parameter/Parameter.h new file mode 100644 index 0000000000000000000000000000000000000000..aec56e55794968995d12eed688988bfd4c284aee --- /dev/null +++ b/ProcessLib/Parameter/Parameter.h @@ -0,0 +1,68 @@ +/** + * \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 PROCESS_LIB_PARAMETER_H_ +#define PROCESS_LIB_PARAMETER_H_ + +#include <memory> +#include <vector> +#include "SpatialPosition.h" + +namespace BaseLib +{ +class ConfigTree; +} // BaseLib + +namespace MeshLib +{ +class Mesh; +} // MeshLib + +namespace ProcessLib +{ +/// Base class for all parameters, not an interface class. This avoids using of +/// void* when storing parameters and convenient destruction. +/// Its property name helps addressing the right parameter. +struct ParameterBase +{ + virtual ~ParameterBase() = default; + + std::string name; +}; + +/*! A Parameter is a function \f$ (t, x) \mapsto f(t, x) \in T^n \f$. + * + * Where \f$ t \f$ is the time and \f$ x \f$ is the SpatialPosition. + * \f$ n \f$ is the number of components of \f$f\f$'s results, i.e., 1 for a + * scalar parameter and >1 for a vectorial or tensorial parameter. + */ +template <typename T> +struct Parameter : public ParameterBase +{ + virtual ~Parameter() = default; + + //! Returns the number of components this Parameter has at every position and + //! point in time. + virtual unsigned getNumberOfComponents() const = 0; + + //! Returns the parameter value at the given time and position. + virtual std::vector<T> const& getTuple( + double const t, SpatialPosition const& pos) const = 0; +}; + +//! Constructs a new ParameterBase from the given configuration. +//! +//! The \c meshes vector is used to set up parameters from mesh input data. +std::unique_ptr<ParameterBase> createParameter( + BaseLib::ConfigTree const& config, + const std::vector<MeshLib::Mesh*>& meshes); + +} // namespace ProcessLib + +#endif // PROCESS_LIB_PARAMETER_H_ diff --git a/ProcessLib/Parameter/SpatialPosition.h b/ProcessLib/Parameter/SpatialPosition.h new file mode 100644 index 0000000000000000000000000000000000000000..fde9bf2f13ed75205c0d69d0122a394a8e99d6fc --- /dev/null +++ b/ProcessLib/Parameter/SpatialPosition.h @@ -0,0 +1,93 @@ +/** + * \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_SPATIALPOSITION_H +#define PROCESSLIB_SPATIALPOSITION_H + +#include <boost/optional.hpp> +#include "MathLib/TemplatePoint.h" + +namespace ProcessLib +{ +//! Represents a position in space which can be either one of +//! a node, an element, an integration point or a cartesian coordinates triple. +//! +//! The setters of this class make sure that only compatible information can be +//! stored at the same time; e.g., it is not possible to specify an element ID +//! and a node ID at the same time (the setAll() method being an exception to +//! that rule). +class SpatialPosition +{ +public: + boost::optional<std::size_t> getNodeID() const { return _node_id; } + boost::optional<std::size_t> getElementID() const { return _element_id; } + boost::optional<unsigned> getIntegrationPoint() const + { + return _integration_point; + } + boost::optional<MathLib::TemplatePoint<double, 3>> const& getCoordinates() + const + { + return _coordinates; + } + + void setNodeID(std::size_t node_id) + { + clear(); + _node_id = node_id; + } + + void setElementID(std::size_t element_id) + { + clear(); + _element_id = element_id; + } + + void setIntegrationPoint(unsigned integration_point) + { + assert(_element_id); + _integration_point = integration_point; + } + + void setCoordinates(MathLib::TemplatePoint<double, 3> const& coordinates) + { + clear(); + _coordinates = coordinates; + } + + void setAll( + boost::optional<std::size_t> const& node_id, + boost::optional<std::size_t> const& element_id, + boost::optional<unsigned> const& integration_point, + boost::optional<MathLib::TemplatePoint<double, 3>> const& coordinates) + { + _node_id = node_id; + _element_id = element_id; + _integration_point = integration_point; + _coordinates = coordinates; + } + + void clear() + { + _node_id = boost::none; + _element_id = boost::none; + _integration_point = boost::none; + _coordinates = boost::none; + } + +private: + boost::optional<std::size_t> _node_id; + boost::optional<std::size_t> _element_id; + boost::optional<unsigned> _integration_point; + boost::optional<MathLib::TemplatePoint<double, 3>> _coordinates; +}; + +} // ProcessLib + +#endif // PROCESSLIB_SPATIALPOSITION_H diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp index d52b3c892b74b1e533d45fb5484e08012583ada1..857e24c9689f68ffc691b585f9d3d637a6506224 100644 --- a/ProcessLib/Process.cpp +++ b/ProcessLib/Process.cpp @@ -69,18 +69,48 @@ void Process::initialize() _process_variables, *_local_to_global_index_map, _integration_order); } -void Process::setInitialConditions(GlobalVector& x) +void Process::setInitialConditions(double const t, 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(); + + 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(t, pos); + + 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,67 +254,12 @@ 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 = NumLib::computeSparsityPattern(*_local_to_global_index_map, _mesh); } -ProcessVariable& findProcessVariable( - std::vector<ProcessVariable> const& variables, - BaseLib::ConfigTree const& pv_config, std::string const& tag) -{ - // Find process variable name in process config. - //! \ogs_file_special - std::string const name = pv_config.getConfigParameter<std::string>(tag); - - // Find corresponding variable by name. - auto variable = std::find_if( - variables.cbegin(), variables.cend(), - [&name](ProcessVariable const& v) { return v.getName() == name; }); - - if (variable == variables.end()) - { - OGS_FATAL( - "Could not find process variable '%s' in the provided variables " - "list for config tag <%s>.", - name.c_str(), tag.c_str()); - } - DBUG("Found process variable \'%s\' for config tag <%s>.", - variable->getName().c_str(), tag.c_str()); - - // Const cast is needed because of variables argument constness. - return const_cast<ProcessVariable&>(*variable); -} - void Process::preIteration(const unsigned iter, const GlobalVector &x) { // In every new iteration cached values of secondary variables are expired. @@ -300,24 +275,4 @@ NumLib::IterationResult Process::postIteration(const GlobalVector &x) return postIterationConcreteProcess(x); } -std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables( - std::vector<ProcessVariable> const& variables, - BaseLib::ConfigTree const& process_config, - std::initializer_list<std::string> - tag_names) -{ - std::vector<std::reference_wrapper<ProcessVariable>> vars; - vars.reserve(tag_names.size()); - - //! \ogs_file_param{process__process_variables} - auto const pv_conf = process_config.getConfigSubtree("process_variables"); - - for (auto const& tag : tag_names) - { - vars.emplace_back(findProcessVariable(variables, pv_conf, tag)); - } - - return vars; -} - } // namespace ProcessLib diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h index b3ca1c3479d98ccf693dfcda5c20b114531aa5df..1f04c0707da1716d7b78009f9d9ceff607c7b9ad 100644 --- a/ProcessLib/Process.h +++ b/ProcessLib/Process.h @@ -15,9 +15,9 @@ #include "NumLib/ODESolver/TimeDiscretization.h" #include "NumLib/NamedFunctionCaller.h" #include "ProcessLib/BoundaryCondition/BoundaryConditionCollection.h" +#include "ProcessLib/Parameter/Parameter.h" #include "ExtrapolatorData.h" -#include "Parameter.h" #include "ProcessOutput.h" #include "SecondaryVariable.h" #include "CachedSecondaryVariable.h" @@ -71,7 +71,7 @@ public: void initialize(); - void setInitialConditions(GlobalVector& x); + void setInitialConditions(const double t, GlobalVector& x); MathLib::MatrixSpecifications getMatrixSpecifications() const override final; @@ -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(); @@ -189,76 +182,6 @@ private: ExtrapolatorData _extrapolator_data; }; -/// Find process variables in \c variables whose names match the settings under -/// the given \c tag_names in the \c process_config. -/// -/// In the process config a process variable is referenced by a name. For -/// example it will be looking for a variable named "H" in the list of process -/// variables when the tag is "hydraulic_head": -/// \code -/// <process> -/// ... -/// <process_variables> -/// <hydraulic_head>H</hydraulic_head> -/// ... -/// </process_variables> -/// ... -/// </process> -/// \endcode -/// -/// \return a vector of references to the found variable(s). -std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables( - std::vector<ProcessVariable> const& variables, - BaseLib::ConfigTree const& process_config, - std::initializer_list<std::string> - tag_names); - -/// Find a parameter of specific type for a name given in the process -/// configuration under the tag. -/// 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... ParameterArgs> -Parameter<ParameterArgs...>& findParameter( - BaseLib::ConfigTree const& process_config, std::string const& tag, - std::vector<std::unique_ptr<ParameterBase>> const& parameters) -{ - // 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; - }); - - 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()); - } - DBUG("Found parameter \'%s\'.", (*parameter_it)->name.c_str()); - - // Check the type correctness of the found parameter. - auto* const parameter = - dynamic_cast<Parameter<ParameterArgs...>*>(parameter_it->get()); - if (!parameter) { - OGS_FATAL("The read parameter is of incompatible type."); - } - return *parameter; -} - } // namespace ProcessLib #endif // PROCESS_LIB_PROCESS_H_ diff --git a/ProcessLib/ProcessVariable.cpp b/ProcessLib/ProcessVariable.cpp index 618c28575347945b0c02863e5614289a5d456f85..6043cf3ee133c4460fa14fdf7f13b656d3bbb8f8 100644 --- a/ProcessLib/ProcessVariable.cpp +++ b/ProcessLib/ProcessVariable.cpp @@ -10,53 +10,29 @@ #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")), + _initial_condition(findParameter<double>( + //! \ogs_file_param{prj__process_variables__process_variable__initial_condition} + config.getConfigParameter<std::string>("initial_condition"), + parameters, _n_components)) { - DBUG("Constructing process variable %s", this->_name.c_str()); - - // Initial condition - //! \ogs_file_param{prj__process_variables__process_variable__initial_condition} - if (auto ic_config = config.getConfigSubtreeOptional("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."); - } - } - else - { - INFO("No initial condition found."); - } + DBUG("Constructing process variable %s", _name.c_str()); // Boundary conditions //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions} diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h index 463c9d1a271df222df1a495af6e84cbdebfef663..2a29dabc3fb860b7aa825dca4f22eca9d0e1dc8a 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; -} +#include "ProcessLib/Parameter/Parameter.h" 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 + Parameter<double> const& getInitialCondition() const { - return _initial_condition->getValue(node_id, component_id); + return _initial_condition; } // Get or create a property vector for results. @@ -78,7 +61,7 @@ private: std::string const _name; MeshLib::Mesh& _mesh; const int _n_components; - std::unique_ptr<InitialCondition> _initial_condition; + Parameter<double> const& _initial_condition; std::vector<BoundaryConditionConfig> _bc_configs; }; diff --git a/ProcessLib/TES/CreateTESProcess.cpp b/ProcessLib/TES/CreateTESProcess.cpp index f14bb0b3af10f080c732a14e3ce95a6e949d5e22..43d7e8ccd8d9b9ccfac9e9150004f7988387dded 100644 --- a/ProcessLib/TES/CreateTESProcess.cpp +++ b/ProcessLib/TES/CreateTESProcess.cpp @@ -9,6 +9,7 @@ #include "CreateTESProcess.h" #include "ProcessLib/Utils/ParseSecondaryVariables.h" +#include "ProcessLib/Utils/ProcessUtils.h" #include "TESProcess.h" namespace ProcessLib diff --git a/ProcessLib/Utils/ProcessUtils.cpp b/ProcessLib/Utils/ProcessUtils.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8fff09c81d6a349975c57fa40676b65bff7287fc --- /dev/null +++ b/ProcessLib/Utils/ProcessUtils.cpp @@ -0,0 +1,62 @@ +/** + * \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 "ProcessUtils.h" +#include "ProcessLib/ProcessVariable.h" + +namespace ProcessLib +{ +ProcessVariable& findProcessVariable( + std::vector<ProcessVariable> const& variables, + BaseLib::ConfigTree const& pv_config, std::string const& tag) +{ + // Find process variable name in process config. + //! \ogs_file_special + std::string const name = pv_config.getConfigParameter<std::string>(tag); + + // Find corresponding variable by name. + auto variable = std::find_if( + variables.cbegin(), variables.cend(), + [&name](ProcessVariable const& v) { return v.getName() == name; }); + + if (variable == variables.end()) + { + OGS_FATAL( + "Could not find process variable '%s' in the provided variables " + "list for config tag <%s>.", + name.c_str(), tag.c_str()); + } + DBUG("Found process variable \'%s\' for config tag <%s>.", + variable->getName().c_str(), tag.c_str()); + + // Const cast is needed because of variables argument constness. + return const_cast<ProcessVariable&>(*variable); +} + +std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables( + std::vector<ProcessVariable> const& variables, + BaseLib::ConfigTree const& process_config, + std::initializer_list<std::string> + tag_names) +{ + std::vector<std::reference_wrapper<ProcessVariable>> vars; + vars.reserve(tag_names.size()); + + //! \ogs_file_param{process__process_variables} + auto const pv_conf = process_config.getConfigSubtree("process_variables"); + + for (auto const& tag : tag_names) + { + vars.emplace_back(findProcessVariable(variables, pv_conf, tag)); + } + + return vars; +} + +} // ProcessLib diff --git a/ProcessLib/Utils/ProcessUtils.h b/ProcessLib/Utils/ProcessUtils.h new file mode 100644 index 0000000000000000000000000000000000000000..b8ac0096e47dbbeaa23bf122779bd6da1d15e3cb --- /dev/null +++ b/ProcessLib/Utils/ProcessUtils.h @@ -0,0 +1,116 @@ +/** + * \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_UTILS_PROCESSUTILS_H +#define PROCESSLIB_UTILS_PROCESSUTILS_H + +#include <vector> +#include "BaseLib/ConfigTree.h" +#include "BaseLib/Error.h" +#include "ProcessLib/Parameter/Parameter.h" + +namespace ProcessLib +{ +class ProcessVariable; + +/// Find process variables in \c variables whose names match the settings under +/// the given \c tag_names in the \c process_config. +/// +/// In the process config a process variable is referenced by a name. For +/// example it will be looking for a variable named "H" in the list of process +/// variables when the tag is "hydraulic_head": +/// \code +/// <process> +/// ... +/// <process_variables> +/// <hydraulic_head>H</hydraulic_head> +/// ... +/// </process_variables> +/// ... +/// </process> +/// \endcode +/// +/// \return a vector of references to the found variable(s). +std::vector<std::reference_wrapper<ProcessVariable>> findProcessVariables( + std::vector<ProcessVariable> const& variables, + BaseLib::ConfigTree const& process_config, + std::initializer_list<std::string> + tag_names); + +/// Find a parameter of specific type for a given name. +/// +/// \see The documentation of the other findParameter() function. +template <typename ParameterDataType> +Parameter<ParameterDataType>& findParameter( + std::string const& parameter_name, + std::vector<std::unique_ptr<ParameterBase>> const& parameters, + unsigned const num_components) +{ + // Find corresponding parameter by 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.", + parameter_name.c_str()); + } + DBUG("Found parameter `%s'.", (*parameter_it)->name.c_str()); + + // Check the type correctness of the found parameter. + auto* const parameter = + dynamic_cast<Parameter<ParameterDataType>*>(parameter_it->get()); + if (!parameter) { + OGS_FATAL("The read parameter `%s' is of incompatible type.", + 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).", + 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 diff --git a/Tests/Data b/Tests/Data index 89e25cfabbc2cbe2ce2dd5493bf200334c70feae..3a9d79f79733b995e38f11f98fdd5fb8ed5d91cd 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit 89e25cfabbc2cbe2ce2dd5493bf200334c70feae +Subproject commit 3a9d79f79733b995e38f11f98fdd5fb8ed5d91cd