[PL] BC; Cleanup assemblies after switch to params

ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h

@@ -54,28 +54,14 @@ public:
         unsigned const n_integration_points =
             Base::_integration_method.getNumberOfPoints();
 
-        SpatialPosition pos;
-        pos.setElementID(id);
-
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-        FemType fe(
-             static_cast<const typename ShapeFunction::MeshElement&>(Base::_element));
-
-        // Get element nodes for the interpolation from nodes to
-        // integration point.
+        // Get element nodes for the interpolation from nodes to integration
+        // point.
         NodalVectorType parameter_node_values =
             _neumann_bc_parameter.getNodalValuesOnElement(Base::_element, t);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
-            pos.setIntegrationPoint(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(ip_data.N));
-            pos.setCoordinates(coords_ip);
 
             _local_rhs.noalias() += ip_data.N *
                                     parameter_node_values.dot(ip_data.N) *

ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h

@@ -21,21 +21,19 @@ namespace ProcessLib
 {
 namespace NormalTractionBoundaryCondition
 {
-template <typename ShapeMatricesTypeDisplacement, int GlobalDim, int NPoints>
+template <typename ShapeMatricesType>
 struct IntegrationPointData final
 {
     IntegrationPointData(
-        typename ShapeMatricesTypeDisplacement::template VectorType<
-            NPoints * GlobalDim> const& Nu_times_n_,
+        typename ShapeMatricesType::ShapeMatrices::ShapeType const N_,
+        typename ShapeMatricesType::GlobalDimVectorType const n_,
         double const integration_weight_)
-        : Nu_times_n(Nu_times_n_), integration_weight(integration_weight_)
+        : N(N_), n(n_), integration_weight(integration_weight_)
     {
     }
 
-    // Shape matrix (for displacement) times element's normal vector.
-    typename ShapeMatricesTypeDisplacement::template VectorType<
-        NPoints * GlobalDim> const Nu_times_n;
-
+    typename ShapeMatricesType::ShapeMatrices::ShapeType const N;
+    typename ShapeMatricesType::GlobalDimVectorType const n;
     double const integration_weight;
 };
 
@@ -56,11 +54,10 @@ class NormalTractionBoundaryConditionLocalAssembler final
     : public NormalTractionBoundaryConditionLocalAssemblerInterface
 {
 public:
-    using ShapeMatricesTypeDisplacement =
+    using ShapeMatricesType =
         ShapeMatrixPolicyType<ShapeFunctionDisplacement, GlobalDim>;
-    using GlobalDimVectorType =
-        typename ShapeMatrixPolicyType<ShapeFunctionDisplacement,
-                                       GlobalDim>::GlobalDimVectorType;
+    using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
 
     NormalTractionBoundaryConditionLocalAssembler(
         MeshLib::Element const& e,
@@ -68,7 +65,9 @@ public:
         bool const is_axially_symmetric,
         unsigned const integration_order,
         Parameter<double> const& pressure)
-        : _integration_method(integration_order), _pressure(pressure)
+        : _integration_method(integration_order),
+          _pressure(pressure),
+          _element(e)
     {
         _local_rhs.setZero(local_matrix_size);
 
@@ -78,10 +77,9 @@ public:
         _ip_data.reserve(n_integration_points);
 
         auto const shape_matrices_u =
-            initShapeMatrices<ShapeFunctionDisplacement,
-                              ShapeMatricesTypeDisplacement, IntegrationMethod,
-                              GlobalDim>(e, is_axially_symmetric,
-                                         _integration_method);
+            initShapeMatrices<ShapeFunctionDisplacement, ShapeMatricesType,
+                              IntegrationMethod, GlobalDim>(
+                e, is_axially_symmetric, _integration_method);
 
         GlobalDimVectorType element_normal(GlobalDim);
 
@@ -104,24 +102,13 @@ public:
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
-            typename ShapeMatricesTypeDisplacement::template MatrixType<
-                GlobalDim, displacement_size>
-                N_u = ShapeMatricesTypeDisplacement::template MatrixType<
-                    GlobalDim, displacement_size>::Zero(GlobalDim,
-                                                        displacement_size);
-            for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
-            {
-                N_u.template block<1, displacement_size / GlobalDim>(
-                       i, i * displacement_size / GlobalDim)
-                    .noalias() = shape_matrices_u[ip].N;
-            }
 
             double const integration_weight =
                 _integration_method.getWeightedPoint(ip).getWeight() *
                 shape_matrices_u[ip].integralMeasure *
                 shape_matrices_u[ip].detJ;
 
-            _ip_data.emplace_back(N_u.transpose() * element_normal,
+            _ip_data.emplace_back(shape_matrices_u[ip].N, element_normal,
                                   integration_weight);
         }
     }
@@ -137,18 +124,28 @@ public:
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();
 
-        SpatialPosition pos;
-        pos.setElementID(id);
+        NodalVectorType pressure =
+            _pressure.getNodalValuesOnElement(_element, t);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
-            pos.setIntegrationPoint(ip);
-
             auto const& w = _ip_data[ip].integration_weight;
-            auto const& Nu_times_n = _ip_data[ip].Nu_times_n;
+            auto const& N = _ip_data[ip].N;
+            auto const& n = _ip_data[ip].n;
 
-            _local_rhs.noalias() -=
-                Nu_times_n.transpose() * _pressure(t, pos)[0] * w;
+            typename ShapeMatricesType::template MatrixType<GlobalDim,
+                                                            displacement_size>
+                N_u = ShapeMatricesType::template MatrixType<
+                    GlobalDim, displacement_size>::Zero(GlobalDim,
+                                                        displacement_size);
+            for (int i = 0; i < static_cast<int>(GlobalDim); ++i)
+            {
+                N_u.template block<1, displacement_size / GlobalDim>(
+                       i, i * displacement_size / GlobalDim)
+                    .noalias() = N;
+            }
+
+            _local_rhs.noalias() -= n.transpose() * N_u * pressure.dot(N) * w;
         }
 
         auto const indices = NumLib::getIndices(id, dof_table_boundary);
@@ -161,17 +158,16 @@ private:
 
     static const int displacement_size =
         ShapeFunctionDisplacement::NPOINTS * GlobalDim;
-    std::vector<IntegrationPointData<ShapeMatricesTypeDisplacement, GlobalDim,
-                                     ShapeFunctionDisplacement::NPOINTS>,
-                Eigen::aligned_allocator<IntegrationPointData<
-                    ShapeMatricesTypeDisplacement, GlobalDim,
-                    ShapeFunctionDisplacement::NPOINTS>>>
+    std::vector<
+        IntegrationPointData<ShapeMatricesType>,
+        Eigen::aligned_allocator<IntegrationPointData<ShapeMatricesType>>>
         _ip_data;
 
-    typename ShapeMatricesTypeDisplacement::template VectorType<
-        displacement_size>
+    typename ShapeMatricesType::template VectorType<displacement_size>
         _local_rhs;
 
+    MeshLib::Element const& _element;
+
 public:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };

ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h

@@ -41,9 +41,7 @@ public:
         bool is_axially_symmetric,
         unsigned const integration_order,
         PythonBoundaryConditionData const& data)
-        : Base(e, is_axially_symmetric, integration_order),
-          _data(data),
-          _element(e)
+        : Base(e, is_axially_symmetric, integration_order), _data(data)
     {
     }
 
@@ -57,12 +55,12 @@ public:
         using FemType =
             NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
         FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &_element));
+            &this->_element));
 
         unsigned const num_integration_points =
             Base::_integration_method.getNumberOfPoints();
         auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
-        auto const num_nodes = _element.getNumberOfNodes();
+        auto const num_nodes = Base::_element.getNumberOfNodes();
         auto const num_comp_total =
             _data.dof_table_bulk.getNumberOfComponents();
 
@@ -85,7 +83,8 @@ public:
                 for (unsigned element_node_id = 0; element_node_id < num_nodes;
                      ++element_node_id)
                 {
-                    auto const* const node = _element.getNode(element_node_id);
+                    auto const* const node =
+                        Base::_element.getNode(element_node_id);
                     auto const boundary_node_id = node->getID();
                     auto const bulk_node_id =
                         bulk_node_ids_map[boundary_node_id];
@@ -119,12 +118,12 @@ public:
 
         for (unsigned ip = 0; ip < num_integration_points; ip++)
         {
-            auto const& sm = Base::_shape_matrices[ip];
-            auto const coords = fe.interpolateCoordinates(sm.N);
+            auto const& N = Base::_ns_and_weights[ip].N;
+            auto const& w = Base::_ns_and_weights[ip].weight;
+            auto const coords = fe.interpolateCoordinates(N);
             prim_vars =
-                sm.N *
-                primary_variables_mat;  // Assumption: all primary variables
-                                        // have same shape functions.
+                N * primary_variables_mat;  // Assumption: all primary variables
+                                            // have same shape functions.
             auto const flag_flux_dFlux =
                 _data.bc_object->getFlux(t, coords, prim_vars_data);
             if (!_data.bc_object->isOverriddenNatural())
@@ -143,9 +142,7 @@ public:
             auto const flux = std::get<1>(flag_flux_dFlux);
             auto const& dFlux = std::get<2>(flag_flux_dFlux);
 
-            auto const& wp = Base::_integration_method.getWeightedPoint(ip);
-            auto const w = sm.detJ * wp.getWeight() * sm.integralMeasure;
-            local_rhs.noalias() += sm.N * (flux * w);
+            local_rhs.noalias() += N * (flux * w);
 
             if (static_cast<int>(dFlux.size()) != num_comp_total)
             {
@@ -170,7 +167,7 @@ public:
                     // The assignement -= takes into account the sign convention
                     // of 1st-order in time ODE systems in OpenGeoSys.
                     local_Jac.block(top, left, width, height).noalias() -=
-                        sm.N.transpose() * (dFlux[comp] * w) * sm.N;
+                        N.transpose() * (dFlux[comp] * w) * N;
                 }
             }
         }
@@ -193,7 +190,6 @@ public:
 
 private:
     PythonBoundaryConditionData const& _data;
-    MeshLib::Element const& _element;
 };
 
 }  // namespace ProcessLib

ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h

@@ -54,22 +54,22 @@ public:
         unsigned const n_integration_points =
             Base::_integration_method.getNumberOfPoints();
 
-        SpatialPosition pos;
-        pos.setElementID(id);
+        typename Base::NodalVectorType const alpha =
+            _data.alpha.getNodalValuesOnElement(Base::_element, t);
+        typename Base::NodalVectorType const u_0 =
+            _data.u_0.getNodalValuesOnElement(Base::_element, t);
 
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
-            pos.setIntegrationPoint(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];
+            auto const& N = ip_data.N;
+            auto const& w = ip_data.weight;
 
             // 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() += ip_data.N * alpha * ip_data.weight;
-            _local_rhs.noalias() += ip_data.N * alpha * u_0 * ip_data.weight;
+            _local_K.diagonal().noalias() += N * alpha.dot(N) * w;
+            _local_rhs.noalias() += N * alpha.dot(N) * u_0.dot(N) * w;
         }
 
         auto const indices = NumLib::getIndices(id, dof_table_boundary);

ProcessLib/BoundaryCondition/VariableDependentNeumannBoundaryConditionLocalAssembler.h

@@ -74,7 +74,8 @@ public:
         NodalVectorType const coefficient_mixed_variables_node_values =
             _data.coefficient_mixed_variables.getNodalValuesOnElement(
                 Base::_element, t);
-        unsigned const n_integration_points = _local_matrix_size;
+        unsigned const n_integration_points =
+            Base::_integration_method.getNumberOfPoints();
 
         auto const indices_current_variable =
             NumLib::getIndices(mesh_item_id, dof_table_boundary);