diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index 4e6d2d880ff4a64b49f0bcaaa426eb225a3c5006..78e8d58f7dadbb4d1fae61499ae92a553d09491f 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -26,6 +26,7 @@ #include "BaseLib/FileTools.h" #include "GeoLib/GEOObjects.h" +#include "MaterialLib/MPL/mpMedium.h" #include "MathLib/Curve/CreatePiecewiseLinearCurve.h" #include "MeshGeoToolsLib/ConstructMeshesFromGeometries.h" #include "MeshGeoToolsLib/CreateSearchLength.h" @@ -233,6 +234,9 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config, //! \ogs_file_param{prj__process_variables} parseProcessVariables(project_config.getConfigSubtree("process_variables")); + //! \ogs_file_param{prj__media} + parseMedia(project_config.getConfigSubtreeOptional("media")); + //! \ogs_file_param{prj__processes} parseProcesses(project_config.getConfigSubtree("processes"), project_directory, output_directory); @@ -309,6 +313,49 @@ void ProjectData::parseParameters(BaseLib::ConfigTree const& parameters_config) parameter->initialize(_parameters); } +void ProjectData::parseMedia( + boost::optional<BaseLib::ConfigTree> const& media_config) +{ + if (!media_config) + return; + + DBUG("Reading media:"); + + if (_mesh_vec.empty() || _mesh_vec[0] == nullptr) + { + ERR("A mesh is required to define medium materials."); + return; + } + + for (auto const& medium_config : + //! \ogs_file_param{material__media__medium} + media_config->getConfigSubtreeList("medium")) + { + auto const material_id = medium_config.getConfigAttribute<int>("id", 0); + + if (_media.find(material_id) != _media.end()) + { + OGS_FATAL( + "Multiple media were specified for the same " + "material id %d. Keep in mind, that if no material id is " + "specified, it is assumed to be 0 by default.", + material_id); + } + + _media[material_id] = + std::make_unique<MaterialPropertyLib::Medium>(medium_config); + } + + if (_media.empty()) + OGS_FATAL("No entity is found inside <media>."); + + if (_media.rbegin()->first != static_cast<int>(_media.size()) - 1) + OGS_FATAL( + "The ids in the porous media definitions in the project file have " + "to be sequential, starting with the id zero."); + +} + void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config, std::string const& project_directory, std::string const& output_directory) diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h index 8fc9248dd189444dcfa1532d09d7886603512644..7dd6472d28ad6856b32dcd1c2cdbf861009c9e39 100644 --- a/Applications/ApplicationsLib/ProjectData.h +++ b/Applications/ApplicationsLib/ProjectData.h @@ -18,7 +18,9 @@ #include <string> #include "BaseLib/ConfigTree.h" +#include "MaterialLib/MPL/mpMedium.h" #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h" + #include "ProcessLib/Parameter/Parameter.h" #include "ProcessLib/Process.h" #include "ProcessLib/ProcessVariable.h" @@ -89,6 +91,9 @@ private: /// Checks if a parameter has name tag. void parseParameters(BaseLib::ConfigTree const& parameters_config); + /// Parses media configuration and saves them in an object. + void parseMedia(boost::optional<BaseLib::ConfigTree> const& media_config); + /// Parses the processes configuration and creates new processes for each /// process entry passing the corresponding subtree to the process /// constructor. @@ -106,6 +111,7 @@ private: void parseCurves(boost::optional<BaseLib::ConfigTree> const& config); + std::unique_ptr<MaterialPropertyLib::Medium> _medium; std::vector<std::unique_ptr<MeshLib::Mesh>> _mesh_vec; std::map<std::string, std::unique_ptr<ProcessLib::Process>> _processes; std::vector<ProcessLib::ProcessVariable> _process_variables; @@ -113,6 +119,8 @@ private: /// Buffer for each parameter config passed to the process. std::vector<std::unique_ptr<ProcessLib::ParameterBase>> _parameters; + std::map<int, std::unique_ptr<MaterialPropertyLib::Medium>> _media; + /// The time loop used to solve this project's processes. std::unique_ptr<TimeLoop> _time_loop; diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt index 2275711f39a896cee82aa70e7376f37346bf57b6..89f4ed54c2866a3cbf9a431e26c8930778f9a08d 100644 --- a/MaterialLib/CMakeLists.txt +++ b/MaterialLib/CMakeLists.txt @@ -13,6 +13,10 @@ append_source_files(SOURCES Fluid/SpecificHeatCapacity) append_source_files(SOURCES Fluid/ThermalConductivity) append_source_files(SOURCES Fluid/WaterVaporProperties) +append_source_files(SOURCES MPL) +append_source_files(SOURCES MPL/Properties) +append_source_files(SOURCES MPL/Components) + append_source_files(SOURCES PorousMedium) append_source_files(SOURCES PorousMedium/Porosity) append_source_files(SOURCES PorousMedium/Storage) diff --git a/MaterialLib/MPL/Components/cWater.cpp b/MaterialLib/MPL/Components/cWater.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cc36f45473dcf35eb942119a5cbd29cc1c908fda --- /dev/null +++ b/MaterialLib/MPL/Components/cWater.cpp @@ -0,0 +1,23 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "cWater.h" +#include "MaterialLib/MPL/Properties/properties.h" + +namespace MaterialPropertyLib +{ +Water::Water() +{ + _properties[name] = std::make_unique<Constant>("Water"); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Components/cWater.h b/MaterialLib/MPL/Components/cWater.h new file mode 100644 index 0000000000000000000000000000000000000000..91279e448de6f3be72def679b0a38b18682a3326 --- /dev/null +++ b/MaterialLib/MPL/Components/cWater.h @@ -0,0 +1,27 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "MaterialLib/MPL/mpComponent.h" + +namespace MaterialPropertyLib +{ +/// A class for Water derived from Component. +/// +/// This class can holds material constants and default properties of ordinary +/// water. +struct Water final : public Component +{ + Water(); +}; +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Components/components.h b/MaterialLib/MPL/Components/components.h new file mode 100644 index 0000000000000000000000000000000000000000..7a330186731f49f3e0b114bff5725f889ad03d5a --- /dev/null +++ b/MaterialLib/MPL/Components/components.h @@ -0,0 +1,15 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "cWater.h" diff --git a/MaterialLib/MPL/CreateMaterialSpatialDistributionMap.cpp b/MaterialLib/MPL/CreateMaterialSpatialDistributionMap.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bece8854a86dfb027442428de9be1082e7c953e9 --- /dev/null +++ b/MaterialLib/MPL/CreateMaterialSpatialDistributionMap.cpp @@ -0,0 +1,48 @@ +/** + * \file + * \date Nov 28, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "CreateMaterialSpatialDistributionMap.h" +#include "MaterialSpatialDistributionMap.h" +#include "MeshLib/Mesh.h" + +namespace MaterialPropertyLib +{ +std::unique_ptr<MaterialSpatialDistributionMap> +createMaterialSpatialDistributionMap( + std::map<int, std::unique_ptr<Medium>> const& media, + MeshLib::Mesh const& mesh) +{ + auto const material_ids = materialIDs(mesh); + + int const max_material_id = + !material_ids + ? 0 + : *std::max_element(begin(*material_ids), end(*material_ids)); + + if (max_material_id > static_cast<int>(media.size() - 1)) + OGS_FATAL( + "The maximum value of MaterialIDs in mesh is %d. As the " + "given number of porous media definitions in the project " + "file is %d, the maximum value of MaterialIDs in mesh must be %d " + "(index starts with zero).", + max_material_id, media.size(), max_material_id - 1); + + if (max_material_id < static_cast<int>(media.size() - 1)) + WARN( + "There are %d porous medium definitions in the project file but " + "only %d different values in the MaterialIDs vector/data_array in " + "the mesh.", + media.size(), max_material_id - 1); + + return std::make_unique<MaterialSpatialDistributionMap>(media, + material_ids); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h b/MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h new file mode 100644 index 0000000000000000000000000000000000000000..97dc2f44b90cc9dc948860fedf0ffd2885fb6eaf --- /dev/null +++ b/MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h @@ -0,0 +1,32 @@ +/** + * \file + * \date Nov 28, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include <map> +#include <memory> + +namespace MeshLib +{ +class Mesh; +} + +namespace MaterialPropertyLib +{ +class MaterialSpatialDistributionMap; + +class Medium; + +std::unique_ptr<MaterialSpatialDistributionMap> +createMaterialSpatialDistributionMap( + std::map<int, std::unique_ptr<Medium>> const& media, + MeshLib::Mesh const& mesh); +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/MaterialSpatialDistributionMap.cpp b/MaterialLib/MPL/MaterialSpatialDistributionMap.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0a81d547026ada9dca1a04da37f0051a40027116 --- /dev/null +++ b/MaterialLib/MPL/MaterialSpatialDistributionMap.cpp @@ -0,0 +1,24 @@ +/** + * \file + * \date Nov 28, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "MaterialSpatialDistributionMap.h" +#include "MeshLib/Mesh.h" + +namespace MaterialPropertyLib +{ +Medium* MaterialSpatialDistributionMap::getMedium(std::size_t const element_id) +{ + auto const material_id = + _material_ids == nullptr ? 0 : (*_material_ids)[element_id]; + + return _media.at(material_id).get(); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/MaterialSpatialDistributionMap.h b/MaterialLib/MPL/MaterialSpatialDistributionMap.h new file mode 100644 index 0000000000000000000000000000000000000000..b49da6d8bd0eeccb2f328edbb799d938e3545ac6 --- /dev/null +++ b/MaterialLib/MPL/MaterialSpatialDistributionMap.h @@ -0,0 +1,44 @@ +/** + * \file + * \date Nov 28, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include <map> +#include <memory> +#include <vector> + +namespace MeshLib +{ +template <typename PROP_VAL_TYPE> +class PropertyVector; +} // namespace MeshLib + +namespace MaterialPropertyLib +{ +class Medium; + +class MaterialSpatialDistributionMap +{ +private: + std::map<int, std::unique_ptr<Medium>> const& _media; + MeshLib::PropertyVector<int> const* const _material_ids; + +public: + MaterialSpatialDistributionMap( + std::map<int, std::unique_ptr<Medium>> const& media, + MeshLib::PropertyVector<int> const* const material_ids) + : _media(media), _material_ids(material_ids) + { + } + + Medium* getMedium(std::size_t element_id); +}; +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/pConstant.cpp b/MaterialLib/MPL/Properties/pConstant.cpp new file mode 100644 index 0000000000000000000000000000000000000000..603ff60cdbe801c138e444751d26c2bd2ef33131 --- /dev/null +++ b/MaterialLib/MPL/Properties/pConstant.cpp @@ -0,0 +1,25 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "pConstant.h" + +namespace MaterialPropertyLib +{ +Constant::Constant(PropertyDataType const& v) +{ + _value = v; + _dvalue = boost::apply_visitor( + [](auto const& value) -> PropertyDataType { return decltype(value){}; }, + v); +}; +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/pConstant.h b/MaterialLib/MPL/Properties/pConstant.h new file mode 100644 index 0000000000000000000000000000000000000000..17bda74f57df56ddddaebdffa9a58e044211ee01 --- /dev/null +++ b/MaterialLib/MPL/Properties/pConstant.h @@ -0,0 +1,30 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "MaterialLib/MPL/mpProperty.h" + +namespace MaterialPropertyLib +{ +/// The constant property class. This property simply retrieves the stored +/// constant value. It accepts all datatypes defined in PropertyDataType +/// (currently: double, Vector, Tensor, std::string) +class Constant final : public Property +{ +public: + /// This constructor accepts single values of any data type defined in the + /// PropertyDataType definition and sets the protected attribute _value of + /// the base class Property to that value. + explicit Constant(PropertyDataType const& v); +}; +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/pUndefined.cpp b/MaterialLib/MPL/Properties/pUndefined.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d4546ce9d2354742ae7ddabf07e4f9468ae1dbb5 --- /dev/null +++ b/MaterialLib/MPL/Properties/pUndefined.cpp @@ -0,0 +1,42 @@ +/** + * \file + * \author Norbert Grunwald + * \date July 3rd, 2018 + * + * \copyright + * Copyright (c) 2012-2018, 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 "pUndefined.h" +#include <string> +#include "MaterialLib/MPL/mpEnums.h" + +namespace MaterialPropertyLib +{ +Undefined::Undefined(PropertyEnum const& pEnum) +{ + thisPropertyEnum = pEnum; +} + +PropertyDataType Undefined::value() const +{ + std::string property = convertEnumToString[thisPropertyEnum]; + + OGS_FATAL( + "The property \'%s\' (property-enum no. %i in " + "MaterialLib/MPL/mpEnums.h) was requested, but is not defined in the " + "project file.", + property.c_str(), (int)thisPropertyEnum); +} + +PropertyDataType Undefined::value(VariableArray const& /*v*/) +{ + return value(); +} +} // MaterialPropertyLib + + diff --git a/MaterialLib/MPL/Properties/pUndefined.h b/MaterialLib/MPL/Properties/pUndefined.h new file mode 100644 index 0000000000000000000000000000000000000000..d614cf30176dd5011a76949d89f88854b084a975 --- /dev/null +++ b/MaterialLib/MPL/Properties/pUndefined.h @@ -0,0 +1,30 @@ +/** + * \file + * \author Norbert Grunwald + * \date July 3rd, 2018 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "MaterialLib/MPL/mpProperty.h" + +namespace MaterialPropertyLib +{ +/// The constant property class. This property simply retrieves the stored +/// constant value. It accepts all datatypes defined in PropertyDataType. +class Undefined final : public Property +{ +private: + PropertyEnum thisPropertyEnum; +public: + explicit Undefined(PropertyEnum const&); + PropertyDataType value() const override; + PropertyDataType value(VariableArray const&) override; +}; +} // MaterialPropertyLib diff --git a/MaterialLib/MPL/Properties/properties.h b/MaterialLib/MPL/Properties/properties.h new file mode 100644 index 0000000000000000000000000000000000000000..0e919e5aa7ad7ab240fe117498b5d8c59660137a --- /dev/null +++ b/MaterialLib/MPL/Properties/properties.h @@ -0,0 +1,16 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "pConstant.h" +#include "pUndefined.h" diff --git a/MaterialLib/MPL/mpComponent.cpp b/MaterialLib/MPL/mpComponent.cpp new file mode 100644 index 0000000000000000000000000000000000000000..940927db78e1565a6825fa33da3b95b1e6bc6695 --- /dev/null +++ b/MaterialLib/MPL/mpComponent.cpp @@ -0,0 +1,84 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "mpComponent.h" +#include "Components/components.h" +#include "Properties/properties.h" + +namespace MaterialPropertyLib +{ +Component::Component() +{ + // Default properties are set to zero. + createDefaultProperties(); + + // Some properties can be initialized by other default properties: + _properties[name] = std::make_unique<Constant>("no_name"); +} +Component::Component(std::string const& component_name) +{ + // Default properties are set to zero. + createDefaultProperties(); + + // Some properties can be initialized by other default properties: + _properties[name] = std::make_unique<Constant>(component_name); +} + +void Component::createProperties(BaseLib::ConfigTree const& config) +{ + for (auto property_config : config.getConfigSubtreeList("property")) + { + // Parsing the property name: + auto const property_name = + property_config.getConfigParameter<std::string>("name"); + // Create a new property based on the configuration subtree: + auto property = newProperty(property_config, this); + + // Insert the new property at the right position into the components + // private PropertyArray: + _properties[convertStringToProperty(property_name)] = + std::move(property); + } +} + +void Component::createDefaultProperties() +{ + for (std::size_t i = 0; i < number_of_property_enums; ++i) + { + _properties[i] = + std::make_unique<Undefined>(static_cast<PropertyEnum>(i)); + } +} + +Property& Component::property(PropertyEnum const& p) const +{ + return *_properties[p]; +} + +std::unique_ptr<Component> newComponent(std::string const& component_name, + bool& isCustomComponent) +{ + // If a name is given, it must conform with one of the derived component + // names in the following list: + if (boost::iequals(component_name, "water")) + { + return std::make_unique<Water>(); + } + + // If the given name does not appear in the list (which is common), the + // method creates a custom component, where all properties have to be + // specified manually. + isCustomComponent = true; + return std::make_unique<Component>(component_name); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpComponent.h b/MaterialLib/MPL/mpComponent.h new file mode 100644 index 0000000000000000000000000000000000000000..9d315653f108bc0236a335c7ecb3c128a2feac3b --- /dev/null +++ b/MaterialLib/MPL/mpComponent.h @@ -0,0 +1,65 @@ +/** + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "BaseLib/ConfigTree.h" +#include "mpProperty.h" + +namespace MaterialPropertyLib +{ +/// \brief This class defines components (substances). +/// \details The Component class is a base class used for not further specified +/// components. Components are specified by the property 'name'. For specified +/// components we derive special classes from this class (for clarity they are +/// located in the 'components' subfolder). +class Component +{ +protected: + /// The property array of the component. + PropertyArray _properties; + +public: + /// Default constructor of Component. This constructor is used + /// when the component is not specified via the 'name'-tag. + Component(); + + /// Constructor for a custom component + explicit Component(std::string const& component_name); + virtual ~Component() = default; + + /// The method reads the 'properties' tag in the prj-file and creates + /// component properties accordingly. + /// + /// First, a new property iy created based on the specified property type. + /// Then, the property name is evaluated and the property is copied into the + /// properties array. + void createProperties(BaseLib::ConfigTree const& /*config*/); + + /// This method initializes the component property array with constant + /// functions of value zero. + void createDefaultProperties(); + /// A get-function for retrieving a cartain property. + Property& property(PropertyEnum const& /*p*/) const; +}; + +/// Method for creating a new component based on the specified component name. +/// +/// This function creates a new component based on the (optional) component name +/// that is given in the prj-file. +/// +/// The method evaluates the string in the 'name'-object and calls the +/// constructors of the derived component classes (if found) or that of the base +/// class (if no name is specified). +std::unique_ptr<Component> newComponent(std::string const& component_name, + bool& isNewComponent); + +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpEnums.h b/MaterialLib/MPL/mpEnums.h new file mode 100644 index 0000000000000000000000000000000000000000..9f185f4c5bde6564220a5a8ca73dae1ec00d7199 --- /dev/null +++ b/MaterialLib/MPL/mpEnums.h @@ -0,0 +1,311 @@ +/** + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#pragma once + +#include <boost/algorithm/string/predicate.hpp> +#include <boost/variant.hpp> +#include <array> +#include <cstddef> +#include <string> +#include "BaseLib/Error.h" + +namespace MaterialPropertyLib +{ +/// Very simple vector data type for holding +/// a pair of values. +using Pair = std::array<double, 2>; + +/// Very simple vector data type for holding +/// vector components. +using Vector = std::array<double, 3>; + +/// Simple symmetric tensor data type for holding +/// xx, yy, zz, xy, xz, yz tensor components. +using SymmTensor = std::array<double, 9>; + +/// Very simple tensor data type for holding +/// tensor components. +using Tensor = std::array<double, 9>; + +/// Variables is simply a list of all commonly used variables that are used to +/// determine the size of the VariableArray. If the variable of your choice is +/// missing, simply add it somewhere at the list, but above the last entry +/// 'numberOfVariables' +enum Variables : std::size_t +{ + phase_pressure, + capillary_pressure, + gas_density, + liquid_density, + temperature, + liquid_saturation, + u, + numberOfVariables +}; + +/// Data type for primary variables, designed to contain both scalar and vector +/// data. +using VariableType = boost::variant<double, Vector>; + +/// The VariableArray is a std::array of fixed size. Its size is determined by +/// the Variables enumerator list. Data type of that array is defined by the +/// VariableType definition. +using VariableArray = std::array<VariableType, numberOfVariables>; + +/// This method returns a value of type double from the variables array +inline double getScalar(VariableType pv) +{ + return boost::get<double>(pv); +} + +/// PropertyEnum is an enumerator list of all known properties of a substance. +/// This includes all properties on all scales (i.e. component, phase, and +/// medium scales). It is used as an indexer for the PropertyArray of the +/// materials. If a necessary property is not in the list, simply add a new one +/// in alphabetical order (of course, except for the last entry +/// 'number_of_property_enums'). Please note that any of these entries must also +/// appear in below convert functions. +enum PropertyEnum : std::size_t +{ + acentric_factor, + binary_interaction_coefficient, + biot_coefficient, + brooks_corey_exponent, + bulk_modulus, + critical_density, + critical_pressure, + critical_temperature, + compressibility, + density, + drhodT, + effective_stress, + entry_pressure, + fredlund_parameters, + heat_capacity, + longitudinal_dispersivity, + molar_mass, + mole_fraction, + molecular_diffusion, + name, + permeability, + phase_velocity, + porosity, + reference_density, + reference_temperature, + reference_pressure, + relative_permeability, + residual_gas_saturation, + residual_liquid_saturation, + saturation, + specific_heat_capacity, + thermal_conductivity, + thermal_expansivity, + transveral_dispersivity, + viscosity, + number_of_property_enums +}; + +/// This function converts a string (e.g. a string from the configuration-tree) +/// into one of the entries of the PropertyEnum enumerator. To avoid confusion, +/// I suggest that the syntax of the properties in the input-files (i.e. the +/// strings) have to be identical to the syntax of the entries in the +/// enumerator. +inline PropertyEnum convertStringToProperty(std::string const& inString) +{ + if (boost::iequals(inString, "acentric_factor")) + { + return acentric_factor; + } + if (boost::iequals(inString, "binary_interaction_coefficient")) + { + return binary_interaction_coefficient; + } + if (boost::iequals(inString, "biot_coefficient")) + { + return biot_coefficient; + } + if (boost::iequals(inString, "brooks_corey_exponent")) + { + return brooks_corey_exponent; + } + if (boost::iequals(inString, "bulk_modulus")) + { + return bulk_modulus; + } + if (boost::iequals(inString, "critical_density")) + { + return critical_density; + } + if (boost::iequals(inString, "critical_pressure")) + { + return critical_pressure; + } + if (boost::iequals(inString, "critical_temperature")) + { + return critical_temperature; + } + if (boost::iequals(inString, "compressibility")) + { + return compressibility; + } + if (boost::iequals(inString, "density")) + { + return density; + } + if (boost::iequals(inString, "drhodT")) + { + return drhodT; + } + if (boost::iequals(inString, "effective_stress")) + { + return effective_stress; + } + if (boost::iequals(inString, "entry_pressure")) + { + return entry_pressure; + } + if (boost::iequals(inString, "fredlund_parameters")) + { + return fredlund_parameters; + } + if (boost::iequals(inString, "heat_capacity")) + { + return heat_capacity; + } + if (boost::iequals(inString, "longitudinal_dispersivity")) + { + return longitudinal_dispersivity; + } + if (boost::iequals(inString, "molar_mass")) + { + return molar_mass; + } + if (boost::iequals(inString, "mole_fraction")) + { + return mole_fraction; + } + if (boost::iequals(inString, "molecular_diffusion")) + { + return molecular_diffusion; + } + if (boost::iequals(inString, "name")) + { + return name; + } + if (boost::iequals(inString, "permeability")) + { + return permeability; + } + if (boost::iequals(inString, "porosity")) + { + return porosity; + } + if (boost::iequals(inString, "phase_velocity")) + { + return phase_velocity; + } + if (boost::iequals(inString, "reference_density")) + { + return reference_density; + } + if (boost::iequals(inString, "reference_temperature")) + { + return reference_temperature; + } + if (boost::iequals(inString, "reference_pressure")) + { + return reference_pressure; + } + if (boost::iequals(inString, "relative_permeability")) + { + return relative_permeability; + } + if (boost::iequals(inString, "residual_gas_saturation")) + { + return residual_gas_saturation; + } + if (boost::iequals(inString, "residual_liquid_saturation")) + { + return residual_liquid_saturation; + } + if (boost::iequals(inString, "saturation")) + { + return saturation; + } + if (boost::iequals(inString, "specific_heat_capacity")) + { + return specific_heat_capacity; + } + if (boost::iequals(inString, "thermal_conductivity")) + { + return thermal_conductivity; + } + if (boost::iequals(inString, "thermal_expansivity")) + { + return thermal_expansivity; + } + if (boost::iequals(inString, "transveral_dispersivity")) + { + return transveral_dispersivity; + } + if (boost::iequals(inString, "viscosity")) + { + return viscosity; + } + + OGS_FATAL( + "The property name \"%s\" does not correspond to any known " + "property", + inString.c_str()); + + return number_of_property_enums; // to avoid the 'no return' warning +} + +const static std::vector<std::string> convertEnumToString{ + {"acentric_factor"}, + {"binary_interaction_coefficient"}, + {"biot_coefficient"}, + {"brooks_corey_exponent"}, + {"bulk_modulus"}, + {"critical_density"}, + {"critical_pressure"}, + {"critical_temperature"}, + {"compressibility"}, + {"density"}, + {"drhodT"}, + {"effective_stress"}, + {"entry_pressure"}, + {"fredlund_parameters"}, + {"heat_capacity"}, + {"longitudinal_dispersivity"}, + {"molar_mass"}, + {"mole_fraction"}, + {"molecular_diffusion"}, + {"name"}, + {"permeability"}, + {"phase_velocity"}, + {"porosity"}, + {"reference_density"}, + {"reference_temperature"}, + {"reference_pressure"}, + {"relative_permeability"}, + {"residual_gas_saturation"}, + {"residual_liquid_saturation"}, + {"saturation"}, + {"specific_heat_capacity"}, + {"thermal_conductivity"}, + {"thermal_expansivity"}, + {"transveral_dispersivity"}, + {"viscosity"}}; + +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpMedium.cpp b/MaterialLib/MPL/mpMedium.cpp new file mode 100644 index 0000000000000000000000000000000000000000..602929c5b864061360dbd7cbcc8fe813986e6ca5 --- /dev/null +++ b/MaterialLib/MPL/mpMedium.cpp @@ -0,0 +1,149 @@ +/** + * \file + * \author Norbert Grunwald + * \date 07.09.2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "mpMedium.h" + +#include <boost/variant.hpp> +#include <string> + +#include "BaseLib/Algorithm.h" + +#include "Properties/pUndefined.h" +#include "Properties/pConstant.h" +#include "mpComponent.h" +#include "mpPhase.h" + +namespace MaterialPropertyLib +{ +Medium::Medium(BaseLib::ConfigTree const& config) +{ + // Default properties are initialized in the first place, such + // that user-defined properties overwrite those defaults. + createDefaultProperties(); + + // A name of a medium is entirely optional and has no effects + // other than a small gain of clarity in case several media + // should be defined. + auto const medium_name = + config.getConfigParameter<std::string>("name", "no_name"); + + // 'name' is a constant property of the medium + _properties[name] = std::make_unique<Constant>(medium_name); + + // Parsing the phases + // Properties of phases may be not required in all the cases. + auto const phases_config = config.getConfigSubtreeOptional("phases"); + if (phases_config) + createPhases(*phases_config); + + // Parsing medium properties, overwriting the defaults. + auto const properties_config = + config.getConfigSubtreeOptional("properties"); + if (properties_config) + { + createProperties(*properties_config); + } + + if (!phases_config && !properties_config) + OGS_FATAL("Neither tag <phases> nor tag <properties> has been found."); +} +void Medium::createPhases(BaseLib::ConfigTree const& config) +{ + std::set<std::string> phase_names; + + for (auto phase_config : config.getConfigSubtreeList("phase")) + { + auto const phase_name = + phase_config.getConfigParameter<std::string>("name"); + + if (phase_name.empty()) + OGS_FATAL("Phase name is a mandatory field and cannot be empty."); + + auto newPhase = std::make_unique<Phase>(phase_name); + // Parsing the components + auto const components_config = + phase_config.getConfigSubtreeOptional("components"); + if (components_config) + newPhase->createComponents(components_config.get()); + + // Properties of phases are optional + auto const properties_config = + phase_config.getConfigSubtreeOptional("properties"); + if (properties_config) + newPhase->createProperties(properties_config.get()); + + if (!components_config && !properties_config) + OGS_FATAL( + "Neither tag <components> nor tag <properties> has been " + "found."); + + phase_names.insert(phase_name); + + _phases.push_back(std::move(newPhase)); + } + + if (phase_names.size() != _phases.size()) + OGS_FATAL( + "Found duplicates with the same phase name tag inside a medium."); +} +void Medium::createProperties(BaseLib::ConfigTree const& config) +{ + for (auto property_config : config.getConfigSubtreeList("property")) + { + // create a new Property based on configuration tree + auto property = newProperty(property_config, this); + /// parse the name of the property + auto const property_name = + property_config.getConfigParameter<std::string>("name"); + // insert the newly created property at the right place + // into the property array + _properties[convertStringToProperty(property_name)] = + std::move(property); + } +} + +void Medium::createDefaultProperties() +{ + for (std::size_t i = 0; i < number_of_property_enums; ++i) + { + this->_properties[i] = std::make_unique<Undefined>((PropertyEnum)i); + } +} + +Phase& Medium::phase(std::size_t const index) const +{ + return *_phases[index]; +} + +Phase& Medium::phase(std::string const& name) const +{ + return *BaseLib::findElementOrError( + _phases.begin(), _phases.end(), + [&name](std::unique_ptr<MaterialPropertyLib::Phase> const& p) { + return getString(p->property( + MaterialPropertyLib::PropertyEnum::name)) == name; + }, + "Could not find phase name '" + name + "'."); +} + +Property& Medium::property(PropertyEnum const& p) const +{ + return *_properties[p]; +} + +std::size_t Medium::numberOfPhases() const +{ + return _phases.size(); +} + +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpMedium.h b/MaterialLib/MPL/mpMedium.h new file mode 100644 index 0000000000000000000000000000000000000000..3737478a9450580200d7147b49c8951684d70702 --- /dev/null +++ b/MaterialLib/MPL/mpMedium.h @@ -0,0 +1,82 @@ +/** + * \author Norbert Grunwald + * \date 07.09.2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include "BaseLib/ConfigTree.h" +#include "mpPhase.h" +#include "mpProperty.h" + +#include <memory> +#include <vector> + +namespace MaterialPropertyLib +{ +/// This class is for material objects on the Medium scale. +/// +/// A Medium consists of an arbitrarily long vector of phases and an array of +/// properties. +class Medium +{ +private: + /// The vector that holds the phases. + std::vector<std::unique_ptr<Phase>> _phases; + /// The array that holds the medium properties. + PropertyArray _properties; + +public: + /// This constructor parses the "phases" and "properties" subtrees of the + /// config tree and calls create methods for the phase vector and the + /// properties array. Medium properties are optional. If not defined, + /// default properties are assigned. + Medium(BaseLib::ConfigTree const& config); + + /// A method that parses the phase details and stores them in the + /// private _phases member. + /// + /// This method creates the phases of the medium. Unlike a medium, a phase + /// may have a name. However, this is silly at the moment since this name + /// still has no effect (except of some benefits in regard of readability). + /// Phase components are required (a phase consists of at least one + /// component). Phase properties are optional. If not given, default + /// properties are assigned. These default properties average the component + /// properties, weighted by mole fraction. + void createPhases(BaseLib::ConfigTree const& config); + + /// A method that parses the medium property details and stores + /// them in the private _properties member. + /// + /// This method creates the properties of the Medium as defined in the + /// prj-file. Only specified properties overwrite the default properties. + void createProperties(BaseLib::ConfigTree const& config); + + /// A method that creates the default properties of the medium. Currently, + /// these defaults is the volume fraction average. + /// + /// Most properties are fine with the volume fraction average, but + /// special-defaults are allowed as well... + void createDefaultProperties(); + + /// A get-function for a particular phase. The ul argument specifies the + /// index within the _phases vector. + Phase& phase(std::size_t index) const; + /// A get-function for a particular phase by phase name. + Phase& phase(std::string const& phase_name) const; + /// A get-function for a property. The argument refers to the name of the + /// property. + Property& property(PropertyEnum const& p) const; + + /// A simple get-function for retrieving the number of phases the medium + /// consists of. + std::size_t numberOfPhases() const; + +}; +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpPhase.cpp b/MaterialLib/MPL/mpPhase.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6fed18b783845569088d71afdd133736d316000c --- /dev/null +++ b/MaterialLib/MPL/mpPhase.cpp @@ -0,0 +1,145 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "mpPhase.h" +#include "BaseLib/Algorithm.h" +#include "Properties/properties.h" + +#include "mpComponent.h" + +namespace MaterialPropertyLib +{ +Phase::Phase(std::string const& phase_name) +{ + createDefaultProperties(); + + std::array<std::string, 4> allowed_phase_names = { + {"solid", "aqueous liquid", "non-aqueous liquid", "gas"}}; + + if (std::find(allowed_phase_names.begin(), + allowed_phase_names.end(), + phase_name) == allowed_phase_names.end()) + { + ERR("Phase name should be one of:"); + for (auto const name : allowed_phase_names) + { + ERR(name.c_str()); + } + OGS_FATAL("Wrong phase name '%s' given.", phase_name.c_str()); + } + _properties[name] = std::make_unique<Constant>(Constant(phase_name)); +} + +void Phase::createComponents(BaseLib::ConfigTree const& config) +{ + // collect names of components + std::set<std::string> component_names; + + for (auto component_config : config.getConfigSubtreeList("component")) + { + // Parsing the component name + auto const component_name = + component_config.getConfigParameter<std::string>("name"); + + // Check whether a name is given or not + if (component_name.empty()) + OGS_FATAL( + "Component name is a mandatory field and cannot be empty."); + + bool isCustomComponent = false; + auto component = newComponent(component_name, isCustomComponent); + + // Parsing component properties. If a component name is given and this + // component is predefined in the class implementation, properties + // become optional. The default values of properties will be overwritten + // if specified. + if (auto const properties_config = + component_config.getConfigSubtreeOptional("properties")) + { + component->createProperties(properties_config.get()); + } + else + { + // No component properties are provided. If the component is + // not specified, this results in a fatal error, since an + // unspecified component has no properties. + if (isCustomComponent) + { + OGS_FATAL("No Properties defined for unspecified component"); + } + } + + component_names.insert(component_name); + _components.push_back(std::move(component)); + } + + if (component_names.size() != _components.size()) + { + OGS_FATAL( + "Found duplicates with the same component name tag inside a " + "phase."); + } +} + +void Phase::createProperties(BaseLib::ConfigTree const& config) +{ + for (auto property_config : config.getConfigSubtreeList("property")) + { + // create a new Property based on configuration tree + auto property = newProperty(property_config, this); + // parse the name of the property + auto const property_name = + property_config.getConfigParameter<std::string>("name"); + // insert the newly created property at the right place into the + // property array + _properties[convertStringToProperty(property_name)] = + std::move(property); + } +} + +void Phase::createDefaultProperties() +{ + for (std::size_t i = 0; i < number_of_property_enums; ++i) + { + _properties[i] = std::make_unique<Undefined>((PropertyEnum)i); + } + + // After this, other special properties can + // be set as exceptional defaults +} + +Component& Phase::component(const std::size_t& index) const +{ + return *_components[index]; +} + +Component& Phase::component(std::string const& name) const +{ + return *BaseLib::findElementOrError( + _components.begin(), _components.end(), + [&name](std::unique_ptr<Component> const& p) { + return getString(p->property(PropertyEnum::name)) == name; + }, + "Could not find component name '" + name + "'."); +} + +Property& Phase::property(PropertyEnum const& p) const +{ + return *_properties[p]; +} + +std::size_t Phase::numberOfComponents() const +{ + return _components.size(); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpPhase.h b/MaterialLib/MPL/mpPhase.h new file mode 100644 index 0000000000000000000000000000000000000000..3dac8cb0cc790af337a9958b674860371294ce55 --- /dev/null +++ b/MaterialLib/MPL/mpPhase.h @@ -0,0 +1,75 @@ +/** + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include <string> +#include <vector> + +#include "BaseLib/ConfigTree.h" + +#include "mpComponent.h" +#include "mpProperty.h" + +namespace MaterialPropertyLib +{ +/// This class defines material phases. +/// +/// The Phase class consists of a vector of components and an array of +/// properties. +class Phase final +{ +private: + std::vector<std::unique_ptr<Component>> _components; + PropertyArray _properties; + +public: + /// The Phase constructor is called with the optional phase name. + explicit Phase(std::string const& phase_name); + + /// The method creating phase components based on config subtree. + /// + /// Just like a phase, a component can have a name. But, in this case, the + /// name has an important task. If a name is given, a specific component + /// class referring to that name with corresponding physical material + /// properties is created. + /// Assigning a name is optional; If no name is given, a custom component + /// without predefined properties is created. + void createComponents(BaseLib::ConfigTree const& config); + + /// The method creating phase properties based on config subtree. + /// + /// This method creates the properties of the Phase as defined in the + /// prj-file. Only specified properties overwrite the default properties. + void createProperties(BaseLib::ConfigTree const& config); + + /// The method initializing the defaule phase properties. + /// Here, all for all properties listed in the Properties enumerator are + /// initialized by mole average functions of value zero. However, + /// 'special-default' properties are allowed to be set. + void createDefaultProperties(); + + /// A simple get-function for a component. The argument refers to the + /// Index of the component in the components vector. + Component& component(std::size_t const& index) const; + + /// A get-function for a component by component name. + Component& component(std::string const& name) const; + + /// A get-function for a property. The argument refers to the name of the + /// property. + Property& property(PropertyEnum const& p) const; + + /// A get-function for retrieving the number of components in this phase. + std::size_t numberOfComponents() const; +}; + +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpProperty.cpp b/MaterialLib/MPL/mpProperty.cpp new file mode 100644 index 0000000000000000000000000000000000000000..331682698b3a8415cf56cb508e426c474c0a25c2 --- /dev/null +++ b/MaterialLib/MPL/mpProperty.cpp @@ -0,0 +1,172 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 "mpProperty.h" + +#include "Properties/properties.h" + +#include "mpComponent.h" +#include "mpMedium.h" +#include "mpPhase.h" + +#include <sstream> +#include <string> +#include <iostream> +#include <vector> + +namespace MaterialPropertyLib +{ + +PropertyDataType Property::value() const +{ + return _value; +} +/// The default implementation of this method only returns the property value +/// without altering it. +PropertyDataType Property::value(VariableArray const&) +{ + return _value; +} + +/// The default implementation of this method only returns the +/// property value derivative without altering it. +PropertyDataType Property::dvalue(VariableArray const&, Variables const) +{ + return _dvalue; +} + +/// Default implementation: 2nd derivative of any constant property is zero. +PropertyDataType Property::ddvalue(VariableArray const&, Variables const, + Variables const) +{ + return 0.0; +} + +void Property::notImplemented(std::string property, std::string material) +{ + OGS_FATAL("The property \'%s\' is not available on the \'%s\' scale", + property.c_str(), material.c_str()); +} + +std::unique_ptr<Property> newProperty(BaseLib::ConfigTree const& config, + Medium* m) +{ + return selectProperty(config, m); +} + +std::unique_ptr<Property> newProperty(BaseLib::ConfigTree const& config, + Phase* p) +{ + return selectProperty(config, p); +} + +std::unique_ptr<Property> newProperty(BaseLib::ConfigTree const& config, + Component* c) +{ + return selectProperty(config, c); +} + +template <typename MaterialType> +std::unique_ptr<Property> selectProperty(BaseLib::ConfigTree const& config, + MaterialType /*material_type*/) +{ + // Parsing the property type: + auto const property_type = config.getConfigParameter<std::string>("type"); + + // If (and only if) the given property type is 'constant', a + // corresponding value is needed. + if (boost::iequals(property_type, "constant")) + { + // TODO (grunwald): Creating constant properties from prj-file is + // currently restricted to scalar values. The Constant constructor, + // however, can handle any datatype defined by PropertyDataType. This + // could be enhanced in order to define vectors or even tensors as + // constant properties. + + auto const sValue = config.getConfigParameter<std::string>("value"); + + std::stringstream issValue(sValue); + + std::vector<double> values; + double dummy(0); + + while (issValue >> dummy) + values.push_back(dummy); + + switch (values.size()) + { + case 1: + { + // scalar + PropertyDataType property_value = values[0]; + return std::make_unique<Constant>(Constant(property_value)); + } + case 2: + { + // Pair + PropertyDataType property_value = (Pair){values[0], values[1]}; + return std::make_unique<Constant>(Constant(property_value)); + } + case 3: + { + // Vector + PropertyDataType property_value = + (Vector){values[0], values[1], values[2]}; + return std::make_unique<Constant>(Constant(property_value)); + } + case 6: + { + // Symmetric Tensor - xx, yy, zz, xy, xz, yz + PropertyDataType property_value = + (SymmTensor){values[0], values[1], values[2], + values[3], values[4], values[5]}; + return std::make_unique<Constant>(Constant(property_value)); + } + case 9: + { + // Tensor + PropertyDataType property_value = + (Tensor){values[0], values[1], values[2], values[3], values[4], + values[5], values[6], values[7], values[8]}; + return std::make_unique<Constant>(Constant(property_value)); + } + + default: + { + OGS_FATAL ("Given number of components for constant property. /%i", values.size()); + } + } + + PropertyDataType property_value; + return std::make_unique<Constant>(Constant(property_value)); + } + // Properties can be medium, phase, or component properties. + // Some of them require information about the respective material. + // Thus, we pass a pointer to the material that requests the property. + // In this method, this pointer is realized via typename MaterialType, which + // replaces either Medium*, Phase*, or Component*. + // Note that most property constructors (only those that request material + // pointers) must be overloaded for any type of material. + + /* TODO Additional properties go here, for example: + if (boost::iequals(property_type, "BilinearTemperaturePressure")) + { + return createBilinearTemperaturePressure(config, material_type); + } + */ + + // If none of the above property types are found, OGS throws an error. + OGS_FATAL("The specified component property type '%s' was not recognized", + property_type.c_str()); +} +} // namespace MaterialPropertyLib diff --git a/MaterialLib/MPL/mpProperty.h b/MaterialLib/MPL/mpProperty.h new file mode 100644 index 0000000000000000000000000000000000000000..e8be00e6f999d572abc25d0b35530f3c2cfb70e1 --- /dev/null +++ b/MaterialLib/MPL/mpProperty.h @@ -0,0 +1,184 @@ +/** + * \file + * \author Norbert Grunwald + * \date Sep 7, 2017 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ +#pragma once + +#include <array> +#include <boost/variant.hpp> +#include <string> + +#include "BaseLib/ConfigTree.h" + +#include "mpEnums.h" + +namespace MaterialPropertyLib +{ +class Medium; +class Phase; +class Component; + +/// This is a custom data type for arbitrary properties, based on the +/// boost::variant container. It can hold scalars, vectors, tensors, and so +/// forth. +enum PropertyDataTypeName +{ + nScalar, + nPair, + nVector, + nSymmTensor, + nTensor +}; + +using PropertyDataType = + boost::variant<double, Pair, Vector, SymmTensor, Tensor, std::string>; + +/// This class is the base class for any material property of any +/// scale (i.e. components, phases, media, ...). The single value of +/// that Property can hold scalars, vectors, tensors, strings, etc. +class Property +{ +protected: + /// The single value of a property. + PropertyDataType _value; + PropertyDataType _dvalue; + +public: + virtual ~Property() = default; + /// This method is called when a property is used for the wrong kind of + /// material, or if the property is not implemented on this kind of material + /// yet. + void notImplemented(std::string /*property*/, std::string /*material*/); + /// This virtual method simply returns the private _value attribute without + /// changing it. + virtual PropertyDataType value() const; + /// This virtual method will compute the property value based on the primary + /// variables that are passed as arguments. + virtual PropertyDataType value(VariableArray const& /*unused*/); + /// This virtual method will compute the derivative of a property + /// with respect to the given variables pv. + virtual PropertyDataType dvalue(VariableArray const&, Variables const pv); + /// This virtual method will compute the second derivative of a + /// property with respect to the given variables pv1 and pv2. + virtual PropertyDataType ddvalue(VariableArray const&, Variables const pv1, + Variables const pv2); +}; + +/// This data type is based on a std::array. It can hold pointers to objects of +/// class Property or its inheritors. The size of this array is determined by +/// the number of entries of the PropertyEnum enumerator. +using PropertyArray = + std::array<std::unique_ptr<Property>, number_of_property_enums>; + +/// Method to select a property by name and to call a derived property +/// constructor. +/// +/// This method creates a new Property. It uses the information stored inside +/// the configuration tree to select the property that is specified by the user. +/// It further passes a pointer to the material that requests the property. +template <typename MaterialType> +std::unique_ptr<Property> selectProperty(BaseLib::ConfigTree const& config, + MaterialType material_type); +/// This method creates a new medium property. +std::unique_ptr<Property> newProperty(BaseLib::ConfigTree const& config, + Medium* m); +/// This method creates a new phase property. +std::unique_ptr<Property> newProperty(BaseLib::ConfigTree const& config, + Phase* p); +/// This method creates a new component property. +std::unique_ptr<Property> newProperty(BaseLib::ConfigTree const& config, + Component* c); +/// This method returns the 0-based index of the variant data types. Can be +/// enhanced by using enums. +inline std::size_t getType(Property const& p) +{ + return p.value().which(); +} + +/// This method returns a value of type double from the property value +/// attribute. +inline double getScalar(Property& p) +{ + assert((getType(p) == PropertyDataTypeName::nScalar) && + "The requested " + "value type is not of type 'double'"); + + return boost::get<double>(p.value()); +} + +/// This method forces the computation of a value of type double and returns it. +inline double getScalar(Property& p, VariableArray const& v) +{ + return boost::get<double>(p.value(v)); +} + +/// This method forces the computation of a value of type Pair and returns it +inline Pair getPair(Property& p, VariableArray const& v) +{ + return boost::get<Pair>(p.value(v)); +} + +/// This method forces the computation of a value of type Vector and returns it. +inline Vector getVector(Property& p) +{ + assert((getType(p) == PropertyDataTypeName::nVector) && + "The requested value type is not of type 'Vector'"); + return boost::get<Vector>(p.value()); +} + +/// This method forces the computation of a value of type Vector and returns it. +inline Vector getVector(Property& p, VariableArray const& v) +{ + return boost::get<Vector>(p.value(v)); +} + +/// This method forces the computation the derivative of a value, returns a pair +/// of double values. The derivative is computed with respect to variable pv. +inline Pair getPairDerivative(Property& p, VariableArray const& v, + Variables const pv) +{ + return boost::get<Pair>(p.dvalue(v, pv)); +} + +/// This method forces the computation the derivative of a value, returns a +/// double value. The derivative is computed with respect to variable pv. +inline double getScalarDerivative(Property& p, VariableArray const& v, + Variables const pv) +{ + return boost::get<double>(p.dvalue(v, pv)); +} + +/// This method forces the computation the second derivative of a value, returns +/// a double value. The derivative is computed with respect to variables pv1 and +/// pv2. +inline double getScalarDerivative(Property& p, VariableArray const& v, + Variables const pv1, Variables const pv2) +{ + return boost::get<double>(p.ddvalue(v, pv1, pv2)); +} + +/// This method returns a value of type string from the property value +/// attribute. +inline std::string getString(Property const& p) +{ + return boost::get<std::string>(p.value()); +} + +/// This method returns a value of any valid type from the property value +/// attribute. The data type is provided by the first parameter in the argument +/// list. +template <typename T> +T getValue(T const&, Property const& p) +{ + return boost::get<T>(p.value()); +} + +} // namespace MaterialPropertyLib diff --git a/Tests/MaterialLib/TestMPL.cpp b/Tests/MaterialLib/TestMPL.cpp new file mode 100644 index 0000000000000000000000000000000000000000..19af610624592435e9cc15f192886c3d133709f9 --- /dev/null +++ b/Tests/MaterialLib/TestMPL.cpp @@ -0,0 +1,24 @@ +/** + * \author Norbert Grunwald + * \date Oct 22, 2018 + * \brief + * + * \copyright + * Copyright (c) 2012-2018, 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 "TestMPL.h" + +MPL::Medium createTestMaterial(std::string const& xml) +{ + auto const ptree = readXml(xml.c_str()); + BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror, + BaseLib::ConfigTree::onwarning); + auto const& config = conf.getConfigSubtree("medium"); + + return MPL::Medium(config); +} diff --git a/Tests/MaterialLib/TestMPL.h b/Tests/MaterialLib/TestMPL.h new file mode 100644 index 0000000000000000000000000000000000000000..7820de6f3c2f1beb4dee653f52c3fb8073be5bcb --- /dev/null +++ b/Tests/MaterialLib/TestMPL.h @@ -0,0 +1,21 @@ +/** + * \file + * \author Norbert Grunwald + * \date Oct 22, 2018 + * + * \copyright + * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#pragma once + +#include "Tests/TestTools.h" + +#include "MaterialLib/MPL/mpMedium.h" +namespace MPL = MaterialPropertyLib; + +MPL::Medium createTestMaterial(std::string const& xml); diff --git a/Tests/MaterialLib/TestMPLParseMaterial.cpp b/Tests/MaterialLib/TestMPLParseMaterial.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d918c9bed9473729fb7dcc65878c822ee911e554 --- /dev/null +++ b/Tests/MaterialLib/TestMPLParseMaterial.cpp @@ -0,0 +1,297 @@ +/** + * \author Norbert Grunwald + * \date Sep 22, 2017 + * + * \copyright + * Copyright (c) 2012-2018, 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 <gtest/gtest.h> + +#include <sstream> + +#include "TestMPL.h" +#include "Tests/TestTools.h" + +#include "MaterialLib/MPL/mpMedium.h" + +namespace MPL = MaterialPropertyLib; + +//---------------------------------------------------------------------------- +// MPL::Medium createTestMaterial(std::string const& xml) +//{ +// auto const ptree = readXml(xml); +// BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror, +// BaseLib::ConfigTree::onwarning); +// auto const& config = +// conf.getConfigSubtree("media").getConfigSubtree("medium"); +// +// return MPL::Medium(config); +//} +/// A simple structure for a string-based component used to create +/// an XML-Structure. +struct Component +{ + std::vector<std::string> property; + Component() : property(MPL::number_of_property_enums) {} +}; +/// A simple structure for a string-based phase used to create +/// an XML-Structure. +struct Phase +{ + std::vector<Component> component; + std::vector<std::string> property; + Phase(std::size_t componentNumber) + : component(componentNumber), property(MPL::number_of_property_enums) + { + } +}; +/// A simple structure for a string-based medium used to create +/// an XML-Structure +struct Medium +{ + std::vector<Phase> phases; + std::vector<std::string> property; + Medium(std::vector<std::size_t> topology) + : property(MPL::number_of_property_enums) + { + for (auto p : topology) + phases.push_back(Phase(p)); + } +}; + +/// Short method creating a number of blanks used for indentation. +std::string indent(std::size_t ind) +{ + return std::string(ind, ' '); +} + +/// This method creates an XML-snippet of a sigle property. +std::string makeProperty(std::size_t ind, std::size_t index, + std::string property) +{ + std::stringstream sProperty; + + if (property == "") + return ""; + + sProperty << indent(ind) << "<property>\n"; + sProperty << indent(ind + 2) << "<name>" << MPL::convertEnumToString[index] + << "</name>\n"; + + // What follows is a C-style string-to-double-conversion... + // seems so easy compared to c++-standard; Feel free to + // change it to latest standard.. + if (double value = atof(property.c_str())) + { + sProperty << indent(ind + 2) << "<type>constant</type>\n"; + sProperty << indent(ind + 2) << "<value>" << value << "</value>\n"; + } + else + sProperty << indent(ind + 2) << "<type>" << property << "</type>\n"; + + sProperty << indent(ind) << "</property>\n"; + return sProperty.str(); +} + +/// This method creates an XML-snippet of a property section. +std::string makeProperties(std::size_t ind, std::vector<std::string> properties) +{ + std::stringstream sProperties; + std::vector<std::string> temp; + bool empty(true); + + for (std::size_t p = 0; p < properties.size(); ++p) + { + std::string property = makeProperty(ind + 2, p, properties[p]); + if ((property != "") && (p != static_cast<std::size_t>(MPL::name))) + { + temp.push_back(property); + empty = false; + } + // If a property is not given, a default takes its place. Note + // also that 'name' is an exceptional property, since it determines + // specific components (thus it is defined outside of the properties- + // tag). + } + + if (empty) + return ""; + sProperties << indent(ind) << "<properties>\n"; + for (auto str : temp) + sProperties << str; + sProperties << indent(ind) << "</properties>\n"; + return sProperties.str(); +} + +/// This method creates an XML-snippet of a single component. +std::string makeComponent(Component c) +{ + std::string componentProperties = makeProperties(14, c.property); + std::stringstream component; + + component << indent(12) << "<component>\n"; + component << indent(14) << "<name>" << c.property[MPL::name] << "</name>\n"; + component << componentProperties; + component << indent(12) << "</component>\n"; + return component.str(); +} + +/// This method creates an XML-snippet of the phase +/// components. +std::string makeComponents(std::vector<Component> components) +{ + std::stringstream sComponents; + sComponents << indent(10) << "<components>\n"; + for (auto c : components) + { + std::string component = makeComponent(c); + sComponents << component; + } + sComponents << indent(10) << "</components>\n"; + return sComponents.str(); +} + +/// This method creates an XML-snippet of a sigle material +/// phase. +std::string makePhase(Phase p) +{ + std::string phaseComponents = makeComponents(p.component); + std::string phaseProperties = makeProperties(8, p.property); + std::stringstream phase; + + phase << indent(8) << "<phase>\n"; + phase << indent(10) << "<name>" << p.property[MPL::name] << "</name>\n"; + phase << phaseComponents; + phase << phaseProperties; + phase << indent(8) << "</phase>\n"; + return phase.str(); +} + +/// This method creates an XML-snippet of the material +/// phases. +std::string makePhases(std::vector<Phase> phases) +{ + std::stringstream sPhases; + + sPhases << indent(6) << "<phases>\n"; + for (auto p : phases) + { + std::string phase = makePhase(p); + sPhases << phase; + } + sPhases << indent(6) << "</phases>\n"; + return sPhases.str(); +} + +/// This method creates the entire XML-tree structure from the string- +/// based medium specification object. I know, indentation is not +/// necessary for this, but my OCD kicked in... +std::string makeMedium(Medium m) +{ + std::string mediumPhases = makePhases(m.phases); + std::string mediumProperties = makeProperties(6, m.property); + std::stringstream medium; + medium << indent(2) << "<medium>\n"; + if (m.property[MPL::name] != "") + medium << indent(4) << "<name>" << m.property[MPL::name] << "</name>\n"; + medium << mediumPhases; + medium << mediumProperties; + medium << indent(2) << "</medium>\n"; + return medium.str(); +} + +/// A method used to obtain the name of a medium, phase, or component of a +/// material or of a specifier and to store them in two vectors for later +/// comparison. +void getNames(MPL::Property const& observation, std::string expectation, + std::string defaultName, std::vector<std::string>* obs, + std::vector<std::string>* exp) +{ + obs->push_back(MPL::getString(observation)); + if (expectation == "") + exp->push_back(defaultName); + else + exp->push_back(expectation); +} + +// This test represents an invariant test. First, several phases, +// components, and properties are generated (more or less randomly). +// Then, an XML-tree is generated from those information. The XML +// tree is then parsed into a configTree, from which an +// MPL::MaterialObject is generated. +// Last step compares the names (as well as the topology) of the +// Material object with the specified parameters. +TEST(Material, parseMaterials) +{ + // This is the topology of our new material: The size of the + // topology vector determines the number of phases, while each + // vector component refers to the number of components of that + // phase. + // The number of properties is fixed in each case and is + // determined by the size of the PropertyEnum enumerator. + std::vector<std::size_t> const mediumTopology = {1, 1, 1}; + + // make a string from every property enumerator + std::vector<std::string> property = MPL::convertEnumToString; + + Medium medium(mediumTopology); + + // the omnivagant medium: + medium.property[MPL::name] = "luminiferous_aether"; + + medium.phases[0].property[MPL::name] = "solid"; + medium.phases[1].property[MPL::name] = "aqueous liquid"; + medium.phases[2].property[MPL::name] = "gas"; + + medium.phases[0].component[0].property[MPL::thermal_conductivity] = "0.654"; + medium.phases[0].component[0].property[MPL::reference_temperature] = "333"; + medium.phases[0].component[0].property[MPL::reference_density] = "2100.0"; + medium.phases[0].component[0].property[MPL::drhodT] = "-0.4"; + + medium.phases[0].component[0].property[MPL::name] = "VerySolid"; + medium.phases[1].component[0].property[MPL::name] = "Water"; + medium.phases[2].component[0].property[MPL::name] = "SuperFluid"; + medium.phases[2].component[0].property[MPL::thermal_conductivity] = "1"; + medium.property[MPL::permeability] = "1.0e-12"; + + // create an actual MaterialProperty-Medium out of the specifier object + auto m = createTestMaterial(makeMedium(medium)); + + // those two vectors will actually be compared + std::vector<std::string> expected; + std::vector<std::string> observed; + + // get the names of the specifier and that of the medium and store + // them in two vectors. If a name of the medium, one of the phases + // or components is not specified, it is automatically replaced by + // the default (which has to be the same as the default of the MPL + getNames(m.property(MPL::name), medium.property[MPL::name], "no_name", + &observed, &expected); + + // now we roam through all phases and components, finding their names + // and storing them in the two vectors + for (std::size_t p = 0; p < m.numberOfPhases(); ++p) + { + const auto& phase = m.phase(p); + getNames(phase.property(MPL::name), + medium.phases[p].property[MPL::name], "no_name", &observed, + &expected); + + for (std::size_t c = 0; c < phase.numberOfComponents(); ++c) + { + const auto& component = phase.component(c); + getNames(component.property(MPL::name), + medium.phases[p].component[c].property[MPL::name], + "no_name", &observed, &expected); + } + } + + // Now, the two vectors are compared. If there is some derivation, + // we can easily locate the problem. + ASSERT_EQ(expected, observed); +}