Parameter.h

Go to the documentation of this file.

 
#pragma once
 
#include <Eigen/Core>
#include <map>
#include <memory>
#include <optional>
#include <utility>
#include <vector>
 
#include "BaseLib/Error.h"
#include "CoordinateSystem.h"
#include "GeoLib/Raster.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Node.h"
#include "SpatialPosition.h"
 
namespace BaseLib
{
class ConfigTree;
}  // namespace BaseLib
 
namespace MathLib
{
class PiecewiseLinearInterpolation;
}  // namespace MathLib
 
namespace MeshLib
{
class Mesh;
}  // namespace MeshLib
 
namespace ParameterLib
{
struct ParameterBase
{
    explicit ParameterBase(std::string name_,
                           MeshLib::Mesh const* mesh = nullptr)
        : name(std::move(name_)), _mesh(mesh)
    {
    }
 
    virtual ~ParameterBase() = default;
 
    virtual bool isTimeDependent() const = 0;
 
    void setCoordinateSystem(CoordinateSystem const& coordinate_system)
    {
        _coordinate_system = coordinate_system;
    }
 
    virtual void initialize(
        std::vector<std::unique_ptr<ParameterBase>> const& /*parameters*/)
    {
    }
 
    MeshLib::Mesh const* mesh() const { return _mesh; }
 
    std::string const name;
 
protected:
    std::vector<double> rotateWithCoordinateSystem(
        std::vector<double> const& values, SpatialPosition const& pos) const
    {
        assert(!!_coordinate_system);  // It is checked before calling this
                                       // function.
 
        // Don't rotate isotropic/scalar values.
        if (values.size() == 1)
        {
            return {values[0]};
        }
        if (values.size() == 2)
        {
            auto const result =
                _coordinate_system->rotateDiagonalTensor<2>(values, pos);
            return {result(0, 0), result(0, 1), result(1, 0), result(1, 1)};
        }
        if (values.size() == 3)
        {
            auto const result =
                _coordinate_system->rotateDiagonalTensor<3>(values, pos);
            return {
                result(0, 0), result(0, 1), result(0, 2),
                result(1, 0), result(1, 1), result(1, 2),
                result(2, 0), result(2, 1), result(2, 2),
            };
        }
        if (values.size() == 4)
        {
            auto const result =
                _coordinate_system->rotateTensor<2>(values, pos);
            return {result(0, 0), result(0, 1), result(1, 0), result(1, 1)};
        }
        if (values.size() == 9)
        {
            auto const result =
                _coordinate_system->rotateTensor<3>(values, pos);
            return {
                result(0, 0), result(0, 1), result(0, 2),
                result(1, 0), result(1, 1), result(1, 2),
                result(2, 0), result(2, 1), result(2, 2),
            };
        }
        OGS_FATAL(
            "Coordinate transformation for a {:d}-component parameter is not "
            "implemented.",
            values.size());
    }
 
protected:
    std::optional<CoordinateSystem> _coordinate_system;
 
    MeshLib::Mesh const* _mesh;
};
 
template <typename T>
struct Parameter : public ParameterBase
{
    using ParameterBase::ParameterBase;
 
    ~Parameter() override = default;
 
    virtual int getNumberOfGlobalComponents() const = 0;
 
    virtual std::vector<T> operator()(double const t,
                                      SpatialPosition const& pos) const = 0;
 
    //
    // The matrix is of the shape NxC, where N is the number of nodes and C is
    // the number of components, such that subsequent multiplication with shape
    // functions matrix from left (which is a row vector) results in a row
    // vector of length C.
    //
    // The default implementation covers all cases, but the derived classes may
    // provide faster implementations.
    virtual Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
    getNodalValuesOnElement(MeshLib::Element const& element,
                            double const t) const
    {
        auto const n_nodes = static_cast<int>(element.getNumberOfNodes());
        auto const n_components = getNumberOfGlobalComponents();
        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> result(n_nodes,
                                                                n_components);
 
        // Column vector of values, copied for each node.
        SpatialPosition x_position;
        auto const nodes = element.getNodes();
        for (int i = 0; i < n_nodes; ++i)
        {
            x_position.setAll(
                nodes[i]->getID(), element.getID(), std::nullopt, *nodes[i]);
            auto const& values = this->operator()(t, x_position);
            auto const row_values =
                Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> const>(
                    values.data(), values.size());
            result.row(i) = row_values;
        }
 
        return result;
    }
};
 
std::unique_ptr<ParameterBase> createParameter(
    BaseLib::ConfigTree const& config,
    const std::vector<std::unique_ptr<MeshLib::Mesh>>& meshes,
    std::vector<GeoLib::NamedRaster> const& named_rasters,
    std::map<std::string,
             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
        curves);
 
std::optional<std::string> isDefinedOnSameMesh(ParameterBase const& parameter,
                                               MeshLib::Mesh const& mesh);
 
}  // namespace ParameterLib