[PL] Unify GenericNeumann and GenericNonuniformN.

These are almost the same and the difference is in the
handling of the shape functions.

ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h

@@ -38,24 +38,58 @@ protected:
     using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
     using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
 
+    struct NAndWeight
+    {
+        NAndWeight(typename ShapeMatricesType::ShapeMatrices::ShapeType&& N_,
+                   double const weight_)
+            : N(std::move(N_)), weight(weight_)
+        {
+        }
+        typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
+        double const weight;
+    };
+
+private:
+    static std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>>
+    initNsAndWeights(MeshLib::Element const& e, bool is_axially_symmetric,
+                     unsigned const integration_order)
+    {
+        IntegrationMethod const integration_method(integration_order);
+        std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>>
+            ns_and_weights;
+        ns_and_weights.reserve(integration_method.getNumberOfPoints());
+
+        auto sms = initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                     IntegrationMethod, GlobalDim>(
+            e, is_axially_symmetric, integration_method);
+        for (unsigned ip = 0; ip < sms.size(); ++ip)
+        {
+            auto& sm = sms[ip];
+            double const w =
+                sm.detJ * sm.integralMeasure *
+                integration_method.getWeightedPoint(ip).getWeight();
+
+            ns_and_weights.emplace_back(std::move(sm.N), w);
+        }
+
+        return ns_and_weights;
+    }
+
 public:
     GenericNaturalBoundaryConditionLocalAssembler(
         MeshLib::Element const& e, bool is_axially_symmetric,
         unsigned const integration_order)
         : _integration_method(integration_order),
-          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                            IntegrationMethod, GlobalDim>(
-              e, is_axially_symmetric, _integration_method)),
+          _ns_and_weights(
+              initNsAndWeights(e, is_axially_symmetric, integration_order)),
           _element(e)
     {
     }
 
 protected:
     IntegrationMethod const _integration_method;
-    std::vector<typename ShapeMatricesType::ShapeMatrices,
-                Eigen::aligned_allocator<
-                    typename ShapeMatricesType::ShapeMatrices>> const
-        _shape_matrices;
+    std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> const
+        _ns_and_weights;
     MeshLib::Element const& _element;
 };
 

ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h

-/**
- * \copyright
- * Copyright (c) 2012-2019, 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 "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
-#include "ProcessLib/Utils/CreateLocalAssemblers.h"
-
-namespace ProcessLib
-{
-template <typename BoundaryConditionData,
-          template <typename, typename, unsigned>
-          class LocalAssemblerImplementation>
-template <typename Data>
-GenericNonuniformNaturalBoundaryCondition<BoundaryConditionData,
-                                          LocalAssemblerImplementation>::
-    GenericNonuniformNaturalBoundaryCondition(
-        unsigned const integration_order, unsigned const shapefunction_order,
-        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
-        int const variable_id, int const component_id,
-        unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data)
-    : _data(std::forward<Data>(data)), _bc_mesh(bc_mesh)
-{
-    static_assert(std::is_same<typename std::decay<BoundaryConditionData>::type,
-                               typename std::decay<Data>::type>::value,
-                  "Type mismatch between declared and passed BC data.");
-
-    // check basic data consistency
-    if (variable_id >=
-            static_cast<int>(dof_table_bulk.getNumberOfVariables()) ||
-        component_id >=
-            dof_table_bulk.getNumberOfVariableComponents(variable_id))
-    {
-        OGS_FATAL(
-            "Variable id or component id too high. Actual values: (%d, %d), "
-            "maximum values: (%d, %d).",
-            variable_id, component_id, dof_table_bulk.getNumberOfVariables(),
-            dof_table_bulk.getNumberOfVariableComponents(variable_id));
-    }
-
-    if (_bc_mesh.getDimension() + 1 != global_dim)
-    {
-        OGS_FATAL(
-            "The dimension of the given boundary mesh (%d) is not by one lower "
-            "than the bulk dimension (%d).",
-            _bc_mesh.getDimension(), global_dim);
-    }
-
-    std::vector<MeshLib::Node*> const& bc_nodes = _bc_mesh.getNodes();
-    DBUG("Found %d nodes for Natural BCs for the variable %d and component %d",
-         bc_nodes.size(), variable_id, component_id);
-
-    MeshLib::MeshSubset bc_mesh_subset(_bc_mesh, bc_nodes);
-
-    // Create local DOF table from the BC mesh subset for the given variable and
-    // component id.
-    _dof_table_boundary.reset(dof_table_bulk.deriveBoundaryConstrainedMap(
-        variable_id, {component_id}, std::move(bc_mesh_subset)));
-
-    createLocalAssemblers<LocalAssemblerImplementation>(
-        global_dim, _bc_mesh.getElements(), *_dof_table_boundary,
-        shapefunction_order, _local_assemblers, _bc_mesh.isAxiallySymmetric(),
-        integration_order, _data);
-}
-
-template <typename BoundaryConditionData,
-          template <typename, typename, unsigned>
-          class LocalAssemblerImplementation>
-void GenericNonuniformNaturalBoundaryCondition<
-    BoundaryConditionData,
-    LocalAssemblerImplementation>::applyNaturalBC(const double t,
-                                                  const GlobalVector& x,
-                                                  GlobalMatrix& K,
-                                                  GlobalVector& b,
-                                                  GlobalMatrix* Jac)
-{
-    GlobalExecutor::executeMemberOnDereferenced(
-        &GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface::
-            assemble,
-        _local_assemblers, *_dof_table_boundary, t, x, K, b, Jac);
-}
-
-}  // namespace ProcessLib

ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition.h

-/**
- * \copyright
- * Copyright (c) 2012-2019, 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 "BoundaryCondition.h"
-#include "MeshLib/MeshSubset.h"
-
-namespace ProcessLib
-{
-class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface;
-
-template <typename BoundaryConditionData,
-          template <typename, typename, unsigned>
-          class LocalAssemblerImplementation>
-class GenericNonuniformNaturalBoundaryCondition final : public BoundaryCondition
-{
-public:
-    /// Create a boundary condition process from given config,
-    /// DOF-table, and a mesh subset for a given variable and its component.
-    /// A local DOF-table, a subset of the given one, is constructed.
-    template <typename Data>
-    GenericNonuniformNaturalBoundaryCondition(
-        unsigned const integration_order, unsigned const shapefunction_order,
-        NumLib::LocalToGlobalIndexMap const& dof_table_bulk,
-        int const variable_id, int const component_id,
-        unsigned const global_dim, MeshLib::Mesh const& bc_mesh, Data&& data);
-
-    /// Calls local assemblers which calculate their contributions to the global
-    /// matrix and the right-hand-side.
-    void applyNaturalBC(const double t, GlobalVector const& x, GlobalMatrix& K,
-                        GlobalVector& b, GlobalMatrix* Jac) override;
-
-private:
-    /// Data used in the assembly of the specific boundary condition.
-    BoundaryConditionData _data;
-
-    /// A lower-dimensional mesh on which the boundary condition is defined.
-    MeshLib::Mesh const& _bc_mesh;
-
-    /// Local dof table, a subset of the global one restricted to the
-    /// participating number of _elements of the boundary condition.
-    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_boundary;
-
-    /// Local assemblers for each element of the boundary mesh.
-    std::vector<std::unique_ptr<
-        GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface>>
-        _local_assemblers;
-};
-
-}  // namespace ProcessLib
-
-#include "GenericNonuniformNaturalBoundaryCondition-impl.h"

ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryConditionLocalAssembler.h

-/**
- * \copyright
- * Copyright (c) 2012-2019, 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 "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-namespace ProcessLib
-{
-class GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface
-{
-public:
-    virtual ~GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface() =
-        default;
-
-    virtual void assemble(
-        std::size_t const id,
-        NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
-        const GlobalVector& x, GlobalMatrix& K, GlobalVector& b,
-        GlobalMatrix* Jac) = 0;
-};
-
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-class GenericNonuniformNaturalBoundaryConditionLocalAssembler
-    : public GenericNonuniformNaturalBoundaryConditionLocalAssemblerInterface
-{
-protected:
-    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
-    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
-    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
-
-    struct NAndWeight
-    {
-        NAndWeight(typename ShapeMatricesType::ShapeMatrices::ShapeType&& N_,
-                   double const weight_)
-            : N(std::move(N_)), weight(weight_)
-        {
-        }
-        typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
-        double const weight;
-    };
-
-private:
-    static std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>>
-    initNsAndWeights(MeshLib::Element const& e, bool is_axially_symmetric,
-                     unsigned const integration_order)
-    {
-        IntegrationMethod const integration_method(integration_order);
-        std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>>
-            ns_and_weights;
-        ns_and_weights.reserve(integration_method.getNumberOfPoints());
-
-        auto sms = initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                     IntegrationMethod, GlobalDim>(
-            e, is_axially_symmetric, integration_method);
-        for (unsigned ip = 0; ip < sms.size(); ++ip)
-        {
-            auto& sm = sms[ip];
-            double const w =
-                sm.detJ * sm.integralMeasure *
-                integration_method.getWeightedPoint(ip).getWeight();
-
-            ns_and_weights.emplace_back(std::move(sm.N), w);
-        }
-
-        return ns_and_weights;
-    }
-
-public:
-    GenericNonuniformNaturalBoundaryConditionLocalAssembler(
-        MeshLib::Element const& e, bool is_axially_symmetric,
-        unsigned const integration_order)
-        : _integration_method(integration_order),
-          _element(e),
-          _ns_and_weights(
-              initNsAndWeights(e, is_axially_symmetric, integration_order))
-    {
-    }
-
-protected:
-    IntegrationMethod const _integration_method;
-    MeshLib::Element const& _element;
-    std::vector<NAndWeight, Eigen::aligned_allocator<NAndWeight>> const
-        _ns_and_weights;
-};
-
-}  // namespace ProcessLib

ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h

@@ -70,15 +70,16 @@ public:
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
-            auto const& sm = Base::_shape_matrices[ip];
+            auto const& ip_data = Base::_ns_and_weights[ip];
             auto const& wp = Base::_integration_method.getWeightedPoint(ip);
 
-            MathLib::TemplatePoint<double, 3> coords_ip(fe.interpolateCoordinates(sm.N));
+            MathLib::TemplatePoint<double, 3> coords_ip(
+                fe.interpolateCoordinates(ip_data.N));
             pos.setCoordinates(coords_ip);
 
-            _local_rhs.noalias() += sm.N * parameter_node_values.dot(sm.N) *
-                                    sm.detJ * wp.getWeight() *
-                                    sm.integralMeasure;
+            _local_rhs.noalias() += ip_data.N *
+                                    parameter_node_values.dot(ip_data.N) *
+                                    ip_data.weight;
         }
 
         auto const indices = NumLib::getIndices(id, dof_table_boundary);

ProcessLib/BoundaryCondition/NonuniformNeumannBoundaryConditionLocalAssembler.h

-/**
- * \copyright
- * Copyright (c) 2012-2019, 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 "MeshLib/PropertyVector.h"
-#include "NumLib/DOF/DOFTableUtil.h"
-#include "NumLib/Function/Interpolation.h"
-#include "ProcessLib/Parameter/MeshNodeParameter.h"
-
-#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
-
-namespace ProcessLib
-{
-struct NonuniformNeumannBoundaryConditionData
-{
-    /* TODO use Parameter */
-    MeshLib::PropertyVector<double> const& values;
-
-    // Used for mapping boundary nodes to bulk nodes.
-    std::size_t bulk_mesh_id;
-    MeshLib::PropertyVector<std::size_t> const& mapping_to_bulk_nodes;
-    NumLib::LocalToGlobalIndexMap const& dof_table_bulk;
-    int const variable_id_bulk;
-    int const component_id_bulk;
-};
-
-template <typename ShapeFunction, typename IntegrationMethod,
-          unsigned GlobalDim>
-class NonuniformNeumannBoundaryConditionLocalAssembler final
-    : public GenericNonuniformNaturalBoundaryConditionLocalAssembler<
-          ShapeFunction, IntegrationMethod, GlobalDim>
-{
-    using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler<
-        ShapeFunction, IntegrationMethod, GlobalDim>;
-    using NodalVectorType = typename Base::NodalVectorType;
-
-public:
-    /// The neumann_bc_value factor is directly integrated into the local
-    /// element matrix.
-    NonuniformNeumannBoundaryConditionLocalAssembler(
-        MeshLib::Element const& e,
-        std::size_t const local_matrix_size,
-        bool const is_axially_symmetric,
-        unsigned const integration_order,
-        NonuniformNeumannBoundaryConditionData const& data)
-        : Base(e, is_axially_symmetric, integration_order),
-          _data(data),
-          _local_rhs(local_matrix_size)
-    {
-    }
-
-    void assemble(std::size_t const id,
-                  NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
-                  double const t, const GlobalVector& /*x*/,
-                  GlobalMatrix& /*K*/, GlobalVector& b,
-                  GlobalMatrix* /*Jac*/) override
-    {
-        _local_rhs.setZero();
-
-        MeshNodeParameter<double> neumann_values{"NeumannValues", _data.values};
-        // Get element nodes for the interpolation from nodes to
-        // integration point.
-        NodalVectorType parameter_node_values =
-            neumann_values.getNodalValuesOnElement(Base::_element, t);
-
-        unsigned const n_integration_points =
-            Base::_integration_method.getNumberOfPoints();
-        for (unsigned ip = 0; ip < n_integration_points; ip++)
-        {
-            auto const& n_and_weight = Base::_ns_and_weights[ip];
-            auto const& N = n_and_weight.N;
-            auto const& w = n_and_weight.weight;
-            _local_rhs.noalias() += N * parameter_node_values.dot(N) * w;
-        }
-
-        auto const indices = NumLib::getIndices(id, dof_table_boundary);
-        b.add(indices, _local_rhs);
-    }
-
-private:
-    NonuniformNeumannBoundaryConditionData const& _data;
-    NodalVectorType _local_rhs;
-
-public:
-    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
-};
-
-}  // namespace ProcessLib

ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h

@@ -60,8 +60,7 @@ public:
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
             pos.setIntegrationPoint(ip);
-            auto const& sm = Base::_shape_matrices[ip];
-            auto const& wp = Base::_integration_method.getWeightedPoint(ip);
+            auto const& ip_data = Base::_ns_and_weights[ip];
 
             double const alpha = _data.alpha(t, pos)[0];
             double const u_0 = _data.u_0(t, pos)[0];
@@ -69,10 +68,8 @@ public:
             // flux = alpha * ( u_0 - u )
             // adding a alpha term to the diagonal of the stiffness matrix
             // and a alpha * u_0 term to the rhs vector
-            _local_K.diagonal().noalias() +=
-                sm.N * alpha * sm.detJ * wp.getWeight() * sm.integralMeasure;
-            _local_rhs.noalias() += sm.N * alpha * u_0 * sm.detJ *
-                                    wp.getWeight() * sm.integralMeasure;
+            _local_K.diagonal().noalias() += ip_data.N * alpha * ip_data.weight;
+            _local_rhs.noalias() += ip_data.N * alpha * u_0 * ip_data.weight;
         }
 
         auto const indices = NumLib::getIndices(id, dof_table_boundary);

ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryCondition.h

@@ -9,14 +9,14 @@
 
 #pragma once
 
-#include "GenericNonuniformNaturalBoundaryCondition.h"
+#include "GenericNaturalBoundaryCondition.h"
 #include "MeshLib/PropertyVector.h"
 #include "VariableDependentNeumannBoundaryConditionLocalAssembler.h"
 
 namespace ProcessLib
 {
 using VariableDependentNeumannBoundaryCondition =
-    GenericNonuniformNaturalBoundaryCondition<
+    GenericNaturalBoundaryCondition<
         VariableDependentNeumannBoundaryConditionData,
         VariableDependentNeumannBoundaryConditionLocalAssembler>;
 

ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h

@@ -14,7 +14,7 @@
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/Parameter/MeshNodeParameter.h"
 
-#include "GenericNonuniformNaturalBoundaryConditionLocalAssembler.h"
+#include "GenericNaturalBoundaryConditionLocalAssembler.h"
 
 namespace ProcessLib
 {
@@ -31,10 +31,10 @@ struct VariableDependentNeumannBoundaryConditionData
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 class VariableDependentNeumannBoundaryConditionLocalAssembler final
-    : public GenericNonuniformNaturalBoundaryConditionLocalAssembler<
+    : public GenericNaturalBoundaryConditionLocalAssembler<
           ShapeFunction, IntegrationMethod, GlobalDim>
 {
-    using Base = GenericNonuniformNaturalBoundaryConditionLocalAssembler<
+    using Base = GenericNaturalBoundaryConditionLocalAssembler<
         ShapeFunction, IntegrationMethod, GlobalDim>;
     using NodalVectorType = typename Base::NodalVectorType;
     using NodalMatrixType = typename Base::NodalMatrixType;