[PL] Impl. of volumetric source terms.

First approach use the original bulk mesh as integration domain.

ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp

+/**
+ * \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 "CreateVolumetricSourceTerm.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/FileTools.h"
+#include "MeshLib/Mesh.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "VolumetricSourceTerm.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
+    BaseLib::ConfigTree const& config,
+    MeshLib::Mesh const& source_term_mesh,
+    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order, unsigned const shapefunction_order,
+    int const variable_id, int const component_id)
+{
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
+    config.checkConfigParameter("type", "Volumetric");
+
+    DBUG("Constructing VolumetricSourceTerm from config.");
+
+    // source term field name
+    auto const& volumetric_source_term_parameter_name =
+        //! \ogs_file_param{prj__process_variables__process_variable_source_terms__volumetric_source_term__parameter}
+        config.getConfigParameter<std::string>("parameter");
+    auto& volumetric_source_term = findParameter<double>(
+        volumetric_source_term_parameter_name, parameters, 1);
+
+    DBUG("Using '%s` as volumetric source term parameter.",
+         volumetric_source_term.name.c_str());
+
+    return std::make_unique<VolumetricSourceTerm>(
+        source_term_mesh, source_term_dof_table, integration_order,
+        shapefunction_order, variable_id, component_id,
+        volumetric_source_term);
+}
+
+}  // namespace ProcessLib

ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h

+/**
+ * \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 <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}
+
+namespace ProcessLib
+{
+struct ParameterBase;
+class SourceTerm;
+
+std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
+    BaseLib::ConfigTree const& config,
+    MeshLib::Mesh const& source_term_mesh,
+    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+    unsigned const integration_order, unsigned const shapefunction_order,
+    int const variable_id, int const component_id);
+
+}   // namespace ProcessLib

ProcessLib/SourceTerms/SourceTerm.h

@@ -17,7 +17,8 @@ namespace ProcessLib
 class SourceTerm
 {
 public:
-    SourceTerm(const NumLib::LocalToGlobalIndexMap& source_term_dof_table)
+    explicit SourceTerm(
+        const NumLib::LocalToGlobalIndexMap& source_term_dof_table)
         : _source_term_dof_table(source_term_dof_table)
     {
     }
@@ -25,6 +26,7 @@ public:
     virtual void integrate(const double t, GlobalVector& b) const = 0;
 
     virtual ~SourceTerm() = default;
+
 protected:
     NumLib::LocalToGlobalIndexMap const& _source_term_dof_table;
 };

ProcessLib/SourceTerms/VolumetricSourceTerm.cpp

+/**
+ * \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 "VolumetricSourceTerm.h"
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+VolumetricSourceTerm::VolumetricSourceTerm(
+    MeshLib::Mesh const& source_term_mesh,
+    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    unsigned const integration_order, unsigned const shapefunction_order,
+    int const variable_id, int const component_id,
+    Parameter<double> const& volumetric_source_term)
+    : SourceTerm(source_term_dof_table),
+      _volumetric_source_term(volumetric_source_term)
+{
+    // check basic data consistency
+    if (variable_id >=
+        static_cast<int>(source_term_dof_table.getNumberOfVariables()))
+    {
+        OGS_FATAL(
+            "Variable id too high. Actual value: %d, maximum value: %d.",
+            variable_id,
+            source_term_dof_table.getNumberOfVariables());
+    }
+    if (component_id >=
+        source_term_dof_table.getNumberOfVariableComponents(variable_id))
+    {
+        OGS_FATAL(
+            "Component id too high. Actual value: %d, maximum value: %d.",
+            component_id,
+            source_term_dof_table.getNumberOfVariableComponents(variable_id));
+    }
+
+    ProcessLib::createLocalAssemblers<VolumetricSourceTermLocalAssembler>(
+        source_term_mesh.getDimension(), source_term_mesh.getElements(),
+        source_term_dof_table, shapefunction_order, _local_assemblers,
+        source_term_mesh.isAxiallySymmetric(), integration_order,
+        _volumetric_source_term);
+}
+
+void VolumetricSourceTerm::integrate(const double t,
+                                     GlobalVector& b) const
+{
+    DBUG("Assemble VolumetricSourceTerm.");
+
+    // Call global assembler for each local assembly item.
+    GlobalExecutor::executeMemberOnDereferenced(
+        &VolumetricSourceTermLocalAssemblerInterface::integrate,
+        _local_assemblers, _source_term_dof_table, t, b);
+}
+
+}   // namespace ProcessLib

ProcessLib/SourceTerms/VolumetricSourceTerm.h

+/**
+ * \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 <memory>
+#include <vector>
+
+#include "SourceTerm.h"
+#include "VolumetricSourceTermFEM.h"
+
+namespace ProcessLib
+{
+class VolumetricSourceTerm final : public SourceTerm
+{
+public:
+    VolumetricSourceTerm(
+        MeshLib::Mesh const& source_term_mesh,
+        NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+        unsigned const integration_order, unsigned const shapefunction_order,
+        int const variable_id, int const component_id,
+        Parameter<double> const& volumetric_source_term);
+
+    void integrate(const double t, GlobalVector& b) const;
+
+private:
+    Parameter<double> const& _volumetric_source_term;
+    std::vector<std::unique_ptr<VolumetricSourceTermLocalAssemblerInterface>>
+        _local_assemblers;
+};
+
+}  // namespace ProcessLib

ProcessLib/SourceTerms/VolumetricSourceTermFEM.h

+/**
+ * \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 <vector>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+class VolumetricSourceTermLocalAssemblerInterface
+{
+public:
+    virtual void integrate(
+        std::size_t const id,
+        NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+        double const t, GlobalVector& b) = 0;
+    virtual ~VolumetricSourceTermLocalAssemblerInterface() = default;
+};
+
+const unsigned NUM_NODAL_DOF = 1;
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+class VolumetricSourceTermLocalAssembler final
+    : public VolumetricSourceTermLocalAssemblerInterface
+{
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+    using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
+        ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+
+    using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
+    using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
+    using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+
+public:
+    VolumetricSourceTermLocalAssembler(
+        MeshLib::Element const& element,
+        std::size_t const local_matrix_size,
+        bool is_axially_symmetric,
+        unsigned const integration_order,
+        Parameter<double> const& volumetric_source_term)
+        : _volumetric_source_term(volumetric_source_term),
+          _integration_method(integration_order),
+          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                            IntegrationMethod, GlobalDim>(
+              element, is_axially_symmetric, _integration_method)),
+          _local_rhs(local_matrix_size)
+    {
+    }
+
+    void integrate(std::size_t const id,
+                   NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+                   double const t, GlobalVector& b) override
+    {
+        _local_rhs.setZero();
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        SpatialPosition pos;
+        pos.setElementID(id);
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            pos.setIntegrationPoint(ip);
+            auto const& sm = _shape_matrices[ip];
+            auto const& wp = _integration_method.getWeightedPoint(ip);
+            auto const st_val = _volumetric_source_term(t, pos)[0];
+
+            _local_rhs.noalias() +=
+                sm.N * st_val * sm.detJ * sm.integralMeasure * wp.getWeight();
+        }
+        auto const indices = NumLib::getIndices(id, source_term_dof_table);
+        b.add(indices, _local_rhs);
+    }
+
+private:
+    Parameter<double> const& _volumetric_source_term;
+
+    IntegrationMethod const _integration_method;
+    std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
+        _shape_matrices;
+    NodalVectorType _local_rhs;
+};
+
+}  // namespace ProcessLib