MaterialLib/MPL/Property.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <Eigen/Core>
#include <array>
#include <string>
#include <typeinfo>
#include <variant>
 
#include "BaseLib/DemangleTypeInfo.h"
#include "BaseLib/Error.h"
#include "ParameterLib/SpatialPosition.h"
#include "PropertyType.h"
#include "VariableType.h"
 
namespace MaterialPropertyLib
{
class Medium;
class Phase;
class Component;
 
using PropertyDataType =
    std::variant<double, Eigen::Matrix<double, 2, 1>,
                 Eigen::Matrix<double, 3, 1>, Eigen::Matrix<double, 2, 2>,
                 Eigen::Matrix<double, 3, 3>, Eigen::Matrix<double, 4, 1>,
                 Eigen::Matrix<double, 6, 1>, Eigen::MatrixXd>;

PropertyDataType fromVector(std::vector<double> const& values);

template <std::size_t N>
PropertyDataType fromArray(std::array<double, N> const& values)
{
    if constexpr (N == 1)
    {
        return values[0];
    }
    else if constexpr (N == 2)
    {
        return Eigen::Vector2d{values[0], values[1]};
    }
    else if constexpr (N == 3)
    {
        return Eigen::Vector3d{values[0], values[1], values[2]};
    }
    else if constexpr (N == 4)
    {
        using M = Eigen::Matrix2d;
        using MRM = Eigen::Matrix<double, 2, 2, Eigen::RowMajor>;
        return M{Eigen::Map<MRM const>{values.data(), 2, 2}};
    }
    else if constexpr (N == 6)
    {
        using M = Eigen::Matrix<double, 6, 1>;
        return M{Eigen::Map<M const>{values.data(), 6}};
    }
    else if constexpr (N == 9)
    {
        using M = Eigen::Matrix3d;
        using MRM = Eigen::Matrix<double, 3, 3, Eigen::RowMajor>;
        return M{Eigen::Map<MRM const>{values.data(), 3, 3}};
    }
    else
    {
        static_assert(false, "Unsupported array size for fromArray.");
    }
}

class Property
{
public:
#ifndef NDEBUG
    virtual ~Property()
    {
        if (property_used)
        {
            DBUG("Property is used: '{:s}'", description());
        }
        else
        {
            WARN("Property is not used: '{:s}'", description());
        }
    }
#else
    virtual ~Property() = default;
#endif

    virtual PropertyDataType initialValue(
        ParameterLib::SpatialPosition const& pos, double const t) const;

    virtual PropertyDataType value() const;
    virtual PropertyDataType value(VariableArray const& variable_array,
                                   VariableArray const& variable_array_prev,
                                   ParameterLib::SpatialPosition const& pos,
                                   double const t, double const dt) const;
    virtual PropertyDataType value(VariableArray const& variable_array,
                                   ParameterLib::SpatialPosition const& pos,
                                   double const t, double const dt) const;
    virtual PropertyDataType dValue(VariableArray const& variable_array,
                                    VariableArray const& variable_array_prev,
                                    Variable const variable,
                                    ParameterLib::SpatialPosition const& pos,
                                    double const t, double const dt) const;
    virtual PropertyDataType dValue(VariableArray const& variable_array,
                                    Variable const variable,
                                    ParameterLib::SpatialPosition const& pos,
                                    double const t, double const dt) const;
    virtual PropertyDataType d2Value(VariableArray const& variable_array,
                                     Variable const variable1,
                                     Variable const variable2,
                                     ParameterLib::SpatialPosition const& pos,
                                     double const t, double const dt) const;

    virtual void setProperties(
        std::vector<std::unique_ptr<Phase>> const& phases);
 
    void setScale(std::variant<Medium*, Phase*, Component*> scale)
    {
        scale_ = scale;
        checkScale();
    };
 
    template <typename T>
    T initialValue(ParameterLib::SpatialPosition const& pos,
                   double const t) const
    {
        try
        {
            return std::get<T>(initialValue(pos, t));
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The initial value of {:s} does not hold requested type '{:s}' "
                "but a {:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_[initialValue(pos, t).index()]);
        }
    }
 
    template <typename T>
    T value() const
    {
        try
        {
#ifndef NDEBUG
            property_used = true;
#endif
            return std::get<T>(value());
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The value of {:s} does not hold requested type '{:s}' but a "
                "{:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_[value().index()]);
        }
    }
 
    template <typename T>
    T value(VariableArray const& variable_array,
            VariableArray const& variable_array_prev,
            ParameterLib::SpatialPosition const& pos, double const t,
            double const dt) const
    {
        try
        {
#ifndef NDEBUG
            property_used = true;
#endif
            return std::get<T>(
                value(variable_array, variable_array_prev, pos, t, dt));
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The value of {:s} is not of the requested type '{:s}' but a "
                "{:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_[value(variable_array,
                                                variable_array_prev, pos, t, dt)
                                              .index()]);
        }
    }
    template <typename T>
    T value(VariableArray const& variable_array,
            ParameterLib::SpatialPosition const& pos, double const t,
            double const dt) const
    {
        try
        {
#ifndef NDEBUG
            property_used = true;
#endif
            return std::get<T>(value(variable_array, pos, t, dt));
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The value of {:s} is not of the requested type '{:s}' but a "
                "{:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_[value(variable_array, pos, t, dt)
                                              .index()]);
        }
    }
 
    template <typename T>
    T dValue(VariableArray const& variable_array,
             VariableArray const& variable_array_prev, Variable const variable,
             ParameterLib::SpatialPosition const& pos, double const t,
             double const dt) const
    {
        try
        {
#ifndef NDEBUG
            property_used = true;
#endif
            return std::get<T>(dValue(variable_array, variable_array_prev,
                                      variable, pos, t, dt));
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The first derivative value of {:s} is not of the requested "
                "type '{:s}' but a {:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_
                    [dValue(variable_array, variable, pos, t, dt).index()]);
        }
    }
    template <typename T>
    T dValue(VariableArray const& variable_array, Variable const variable,
             ParameterLib::SpatialPosition const& pos, double const t,
             double const dt) const
    {
        try
        {
#ifndef NDEBUG
            property_used = true;
#endif
            return std::get<T>(dValue(variable_array, variable, pos, t, dt));
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The first derivative value of {:s} is not of the requested "
                "type '{:s}' but a {:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_
                    [dValue(variable_array, variable, pos, t, dt).index()]);
        }
    }
    template <typename T>
    T d2Value(VariableArray const& variable_array, Variable const& variable1,
              Variable const& variable2,
              ParameterLib::SpatialPosition const& pos, double const t,
              double const dt) const
    {
        try
        {
#ifndef NDEBUG
            property_used = true;
#endif
            return std::get<T>(
                d2Value(variable_array, variable1, variable2, pos, t, dt));
        }
        catch (std::bad_variant_access const&)
        {
            OGS_FATAL(
                "The second derivative value of {:s} is not of the requested "
                "type '{:s}' but a {:s}.",
                description(),
                BaseLib::typeToString<T>(),
                property_data_type_names_[d2Value(variable_array, variable1,
                                                  variable2, pos, t, dt)
                                              .index()]);
        }
    }
 
protected:
    std::string name_;
    PropertyDataType value_;
    PropertyDataType dvalue_;
    std::variant<Medium*, Phase*, Component*> scale_;
 
private:
    virtual void checkScale() const
    {
        // Empty check for properties which can be defined on every scale,
        // medium, phase or component
    }
    std::string description() const;
#ifndef NDEBUG
    mutable bool property_used = false;
#endif
 
private:
    static constexpr std::array property_data_type_names_ = {
        "scalar",           "2-vector",           "3-vector",
        "2x2-matrix",       "3x3-matrix",         "2D-Kelvin vector",
        "3D-Kelvin vector", "dynamic matrix type"};
    static_assert(property_data_type_names_.size() ==
                      std::variant_size_v<PropertyDataType>,
                  "The array of property data type names has different size "
                  "than the PropertyDataType variant type.");
};
 
inline void overwriteExistingProperties(
    PropertyArray& properties,
    PropertyArray& new_properties,
    std::variant<Medium*, Phase*, Component*>
        scale_pointer)
{
    for (std::size_t i = 0; i < properties.size(); ++i)
    {
        if (new_properties[i] != nullptr)
        {
            properties[i] = std::move(new_properties[i]);
            properties[i]->setScale(scale_pointer);
        }
    }
}
 
inline void updatePropertiesForAllPhases(
    PropertyArray& properties,
    std::vector<std::unique_ptr<Phase>> const& phases)
{
    for (auto& p : properties)
    {
        if (p != nullptr)
        {
            p->setProperties(phases);
        }
    }
}
 
}  // namespace MaterialPropertyLib