[PL] added RandomFieldMeshElementParameter

Documentation/ProjectFile/prj/parameters/parameter/RandomFieldMeshElementParameter/c_RandomFieldMeshElementParameter.md

+Defines a parameter with random values in each mesh element with a uniform distribution in a given range. A given seed for initializing the random number generator allows for repeatability of the random field.

Documentation/ProjectFile/prj/parameters/parameter/RandomFieldMeshElementParameter/t_field_name.md

+The name of the mesh property in which the parameter values are written.

Documentation/ProjectFile/prj/parameters/parameter/RandomFieldMeshElementParameter/t_range.md

+Defines the minimum and maximum value of the uniform distribution, specified by two floats.

Documentation/ProjectFile/prj/parameters/parameter/RandomFieldMeshElementParameter/t_seed.md

+The seed to initialize the random number generator.

ParameterLib/Parameter.cpp

@@ -18,6 +18,7 @@
 #include "GroupBasedParameter.h"
 #include "MeshElementParameter.h"
 #include "MeshNodeParameter.h"
+#include "RandomFieldMeshElementParameter.h"
 #include "TimeDependentHeterogeneousParameter.h"
 
 namespace ParameterLib
@@ -76,6 +77,15 @@ std::unique_ptr<ParameterBase> createParameter(
         INFO("MeshNodeParameter: {:s}", name);
         return createMeshNodeParameter(name, config, mesh);
     }
+    if (type == "RandomFieldMeshElement")
+    {
+        auto& mesh_var = *BaseLib::findElementOrError(
+            begin(meshes), end(meshes),
+            [&mesh_name](auto const& m) { return m->getName() == mesh_name; },
+            "Expected to find a mesh named " + mesh_name + ".");
+        INFO("RandomFieldMeshElement: {:s}", name);
+        return createRandomFieldMeshElementParameter(name, config, mesh_var);
+    }
     if (type == "TimeDependentHeterogeneousParameter")
     {
         INFO("TimeDependentHeterogeneousParameter: {:s}", name);

ParameterLib/RandomFieldMeshElementParameter.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "RandomFieldMeshElementParameter.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "MeshLib/Mesh.h"
+#include <random>
+#include <functional>
+
+namespace ParameterLib
+{
+std::unique_ptr<ParameterBase> createRandomFieldMeshElementParameter(
+    std::string const& name, BaseLib::ConfigTree const& config,
+    MeshLib::Mesh& mesh)
+{
+    //! \ogs_file_param{prj__parameters__parameter__type}
+    config.checkConfigParameter("type", "RandomFieldMeshElement");
+    auto const field_name =
+        //! \ogs_file_param{prj__parameters__parameter__RandomFieldMeshElement__field_name}
+        config.getConfigParameter<std::string>("field_name");
+    auto const range =
+        //! \ogs_file_param{prj__parameters__parameter__RandomFieldMeshElement__range}
+        config.getConfigParameter<std::vector<double>>("range");
+    if (range.size() != 2)
+    {
+        OGS_FATAL(
+            "The range needs to have two components, but {:d} were given.",
+            range.size());
+    }
+    auto const seed =
+        //! \ogs_file_param{prj__parameters__parameter__RandomFieldMeshElement__seed}
+        config.getConfigParameter<int>("seed");
+    DBUG("Generating field {:s} with range {:g} to {:g} and seed {:d}.",
+         field_name, range[0], range[1], seed);
+
+    std::vector<double> values(mesh.getElements().size());
+
+    std::mt19937 generator(seed);
+    std::uniform_real_distribution<> distr(range[0], range[1]);
+    auto gen = [&distr, &generator]() { return distr(generator); };
+    generate(begin(values), end(values), gen);
+
+    MeshLib::addPropertyToMesh(mesh, field_name, MeshLib::MeshItemType::Cell, 1,
+                               values);
+
+    auto const& property =
+        mesh.getProperties().getPropertyVector<double>(field_name);
+
+    if (property->getMeshItemType() != MeshLib::MeshItemType::Cell)
+    {
+        OGS_FATAL("The mesh property `{:s}' is not an element property.",
+                  field_name);
+    }
+
+    return std::make_unique<RandomFieldMeshElementParameter<double>>(name, mesh,
+                                                                     *property);
+}
+
+}  // namespace ParameterLib

ParameterLib/RandomFieldMeshElementParameter.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "Parameter.h"
+
+namespace MeshLib
+{
+template <typename T>
+class PropertyVector;
+}  // namespace MeshLib
+
+namespace ParameterLib
+{
+/// A parameter represented by a mesh property vector.
+template <typename T>
+struct RandomFieldMeshElementParameter final : public Parameter<T>
+{
+    RandomFieldMeshElementParameter(std::string const& name_,
+                                    MeshLib::Mesh& mesh,
+                                    MeshLib::PropertyVector<T> const& property)
+        : Parameter<T>(name_, &mesh), _property(property)
+    {
+    }
+
+    bool isTimeDependent() const override { return false; }
+
+    int getNumberOfGlobalComponents() const override
+    {
+        return _property.getNumberOfGlobalComponents();
+    }
+
+    std::vector<T> operator()(double const /*t*/,
+                              SpatialPosition const& pos) const override
+    {
+        auto const e = pos.getElementID();
+        if (!e)
+        {
+            OGS_FATAL(
+                "Trying to access a RandomFieldMeshElementParameter but "
+                "the element id is not specified.");
+        }
+        auto const num_comp = _property.getNumberOfGlobalComponents();
+        std::vector<T> cache(num_comp);
+        for (int c = 0; c < num_comp; ++c)
+        {
+            cache[c] = _property.getComponent(*e, c);
+        }
+
+        if (!this->_coordinate_system)
+        {
+            return cache;
+        }
+
+        return this->rotateWithCoordinateSystem(cache, pos);
+    }
+
+    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> getNodalValuesOnElement(
+        MeshLib::Element const& element, double const t) const override
+    {
+        auto const n_nodes = element.getNumberOfNodes();
+        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> result(
+            n_nodes, getNumberOfGlobalComponents());
+
+        // Column vector of values, copied for each node.
+        SpatialPosition x_position;
+        x_position.setElementID(element.getID());
+        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());
+        for (unsigned i = 0; i < n_nodes; ++i)
+        {
+            result.row(i) = row_values;
+        }
+        return result;
+    }
+
+private:
+    MeshLib::PropertyVector<T> const& _property;
+};
+
+std::unique_ptr<ParameterBase> createRandomFieldMeshElementParameter(
+    std::string const& name, BaseLib::ConfigTree const& config,
+    MeshLib::Mesh& mesh);
+
+}  // namespace ParameterLib