[PL/Stokes] add local assembler.

ProcessLib/StokesFlow/CreateLocalAssemblers.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "LocalDataInitializer.h"
+
+#include "BaseLib/Logging.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+namespace ProcessLib
+{
+namespace StokesFlow
+{
+namespace detail
+{
+template <int GlobalDim,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    const unsigned shapefunction_order,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    // Shape matrices initializer
+    using LocalDataInitializer =
+        LocalDataInitializer<LocalAssemblerInterface,
+                             LocalAssemblerImplementation, GlobalDim,
+                             ExtraCtorArgs...>;
+
+    DBUG("Create local assemblers.");
+    // Populate the vector of local assemblers.
+    local_assemblers.resize(mesh_elements.size());
+
+    LocalDataInitializer initializer(dof_table, shapefunction_order);
+
+    DBUG("Calling local assembler builder for all mesh elements.");
+    GlobalExecutor::transformDereferenced(
+        initializer, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+}  // namespace detail
+
+/*! Creates local assemblers for each element of the given \c mesh.
+ *
+ * \tparam LocalAssemblerImplementation the individual local assembler type
+ * \tparam LocalAssemblerInterface the general local assembler interface
+ * \tparam ExtraCtorArgs types of additional constructor arguments.
+ *         Those arguments will be passed to the constructor of
+ *         \c LocalAssemblerImplementation.
+ *
+ * The first two template parameters cannot be deduced from the arguments.
+ * Therefore they always have to be provided manually.
+ */
+template <int GlobalDim,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblers(
+    const unsigned /*dimension*/,
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    const unsigned shapefunction_order,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    DBUG("Create local assemblers.");
+
+    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
+        dof_table, shapefunction_order, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+}  // namespace StokesFlow
+}  // namespace ProcessLib

ProcessLib/StokesFlow/LocalDataInitializer.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <functional>
+#include <memory>
+#include <type_traits>
+#include <typeindex>
+#include <typeinfo>
+#include <unordered_map>
+
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/FiniteElement/LowerDimShapeTable.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+
+#ifndef OGS_MAX_ELEMENT_DIM
+static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
+#endif
+
+#ifndef OGS_MAX_ELEMENT_ORDER
+static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
+#endif
+
+// The following macros decide which element types will be compiled, i.e.
+// which element types will be available for use in simulations.
+
+#ifdef OGS_ENABLE_ELEMENT_SIMPLEX
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 1u
+#else
+#define ENABLED_ELEMENT_TYPE_SIMPLEX 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_CUBOID
+#define ENABLED_ELEMENT_TYPE_CUBOID 1u << 1
+#else
+#define ENABLED_ELEMENT_TYPE_CUBOID 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PRISM
+#define ENABLED_ELEMENT_TYPE_PRISM 1u << 2
+#else
+#define ENABLED_ELEMENT_TYPE_PRISM 0u
+#endif
+
+#ifdef OGS_ENABLE_ELEMENT_PYRAMID
+#define ENABLED_ELEMENT_TYPE_PYRAMID 1u << 3
+#else
+#define ENABLED_ELEMENT_TYPE_PYRAMID 0u
+#endif
+
+// Dependent element types.
+// Faces of tets, pyramids and prisms are triangles
+#define ENABLED_ELEMENT_TYPE_TRI                                       \
+    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+     (ENABLED_ELEMENT_TYPE_PRISM))
+// Faces of hexes, pyramids and prisms are quads
+#define ENABLED_ELEMENT_TYPE_QUAD                                     \
+    ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
+     (ENABLED_ELEMENT_TYPE_PRISM))
+
+// All enabled element types
+#define OGS_ENABLED_ELEMENTS                                          \
+    ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
+     (ENABLED_ELEMENT_TYPE_PYRAMID) | (ENABLED_ELEMENT_TYPE_PRISM))
+
+// Include only what is needed (Well, the conditions are not sharp).
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0
+#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
+#endif
+
+namespace ProcessLib::StokesFlow
+{
+/// The LocalDataInitializer is a functor creating a local assembler data with
+/// corresponding to the mesh element type shape functions and calling
+/// initialization of the new local assembler data.
+/// For example for MeshLib::Line a local assembler data with template argument
+/// NumLib::ShapeLine2 is created.
+///
+/// \attention This is modified version of the ProcessLib::LocalDataInitializer
+/// class which does not include line or point elements. For the shape functions
+/// of order 2 (used for fluid velocity in StokesFlow) a shape function of order
+/// 1 will be used for the pressure.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, typename, int> class LocalAssemblerData,
+          unsigned GlobalDim, typename... ConstructorArgs>
+class LocalDataInitializer final
+{
+public:
+    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+
+    LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table,
+                         const unsigned shapefunction_order)
+        : _dof_table(dof_table)
+    {
+        if (shapefunction_order != 2)
+        {
+            OGS_FATAL(
+                "The given shape function order {:d} is not supported. The "
+                "shape function order shall be specified to 2.",
+                shapefunction_order);
+        }
+
+        // /// Quads and Hexahedra ///////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+            _builder[std::type_index(typeid(MeshLib::Quad8))] =
+                makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
+            _builder[std::type_index(typeid(MeshLib::Quad9))] =
+                makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+            _builder[std::type_index(typeid(MeshLib::Hex20))] =
+                makeLocalAssemblerBuilder<NumLib::ShapeHex20>();
+#endif
+
+            // /// Simplices ////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
+            _builder[std::type_index(typeid(MeshLib::Tri6))] =
+                makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
+#endif
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+            _builder[std::type_index(typeid(MeshLib::Tet10))] =
+                makeLocalAssemblerBuilder<NumLib::ShapeTet10>();
+#endif
+
+            // /// Prisms ////////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+            _builder[std::type_index(typeid(MeshLib::Prism15))] =
+                makeLocalAssemblerBuilder<NumLib::ShapePrism15>();
+#endif
+
+            // /// Pyramids //////////////////////////////////////////////////
+
+#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PYRAMID) != 0 && \
+    OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
+            _builder[std::type_index(typeid(MeshLib::Pyramid13))] =
+                makeLocalAssemblerBuilder<NumLib::ShapePyra13>();
+#endif
+    }
+
+    /// Returns data pointer to the newly created local assembler data.
+    ///
+    /// \attention
+    /// The index \c id is not necessarily the mesh item's id. Especially when
+    /// having multiple meshes it will differ from the latter.
+    LADataIntfPtr operator()(std::size_t const id,
+                             MeshLib::Element const& mesh_item,
+                             ConstructorArgs&&... args) const
+    {
+        auto const type_idx = std::type_index(typeid(mesh_item));
+        auto const it = _builder.find(type_idx);
+
+        if (it != _builder.end())
+        {
+            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
+            return it->second(mesh_item, num_local_dof,
+                              std::forward<ConstructorArgs>(args)...);
+        }
+        OGS_FATAL(
+            "You are trying to build a local assembler for an unknown mesh "
+            "element type ({:s})."
+            " Maybe you have disabled this mesh element type in your build "
+            "configuration.",
+            type_idx.name());
+    }
+
+private:
+    using LADataBuilder =
+        std::function<LADataIntfPtr(MeshLib::Element const& e,
+                                    std::size_t const local_matrix_size,
+                                    ConstructorArgs&&...)>;
+
+    template <typename ShapeFunction>
+    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+        typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+    template <typename ShapeFunctionDisplacement,
+              typename ShapeFunctionPressure>
+    using LAData =
+        LocalAssemblerData<ShapeFunctionDisplacement, ShapeFunctionPressure,
+                           IntegrationMethod<ShapeFunctionDisplacement>,
+                           GlobalDim>;
+
+    /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
+    /// depending whether the global dimension is less than the shape function's
+    /// dimension or not.
+    template <typename ShapeFunction>
+    static LADataBuilder makeLocalAssemblerBuilder()
+    {
+        return makeLocalAssemblerBuilder<ShapeFunction>(
+            static_cast<std::integral_constant<
+                bool, (GlobalDim >= ShapeFunction::DIM)>*>(nullptr));
+    }
+
+    /// Mapping of element types to local assembler constructors.
+    std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+    // local assembler builder implementations.
+private:
+    /// Generates a function that creates a new LocalAssembler of type
+    /// LAData<ShapeFunctionDisplacement, ShapeFunctionPressure>. Only functions
+    /// with shape function's dimension less or equal to the global dimension
+    /// are instantiated, e.g.  following combinations of shape functions and
+    /// global dimensions: (Line2, 1), (Line2, 2), (Line2, 3), (Hex20, 3) but
+    /// not (Hex20, 2) or (Hex20, 1).
+    template <typename ShapeFunction>
+    static LADataBuilder makeLocalAssemblerBuilder(std::true_type* /*unused*/)
+    {
+        assert(ShapeFunction::ORDER == 2);
+
+        using LowerOrderShapeFunction =
+            typename NumLib::LowerDim<ShapeFunction>::type;
+        return [](MeshLib::Element const& e,
+                  std::size_t const local_matrix_size,
+                  ConstructorArgs&&... args) {
+            return LADataIntfPtr{
+                new LAData<ShapeFunction, LowerOrderShapeFunction>{
+                    e, local_matrix_size,
+                    std::forward<ConstructorArgs>(args)...}};
+        };
+    }
+
+    /// Returns nullptr for shape functions whose dimensions are less than the
+    /// global dimension.
+    template <typename ShapeFunction>
+    static LADataBuilder makeLocalAssemblerBuilder(std::false_type* /*unused*/)
+    {
+        return nullptr;
+    }
+};
+
+}  // namespace ProcessLib::StokesFlow
+
+#undef ENABLED_ELEMENT_TYPE_SIMPLEX
+#undef ENABLED_ELEMENT_TYPE_CUBOID
+#undef ENABLED_ELEMENT_TYPE_PYRAMID
+#undef ENABLED_ELEMENT_TYPE_PRISM
+#undef ENABLED_ELEMENT_TYPE_TRI
+#undef ENABLED_ELEMENT_TYPE_QUAD
+#undef OGS_ENABLED_ELEMENTS

ProcessLib/StokesFlow/StokesFlowFEM.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <Eigen/Dense>
+#include <vector>
+
+#include "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
+#include "StokesFlowProcessData.h"
+
+#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Property.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "NumLib/Function/Interpolation.h"
+
+namespace ProcessLib::StokesFlow
+{
+template <typename ShapeFunctionLiquidVelocity, typename ShapeFunctionPressure,
+          typename IntegrationMethod, int GlobalDim>
+class LocalAssemblerData : public StokesFlowLocalAssemblerInterface
+{
+    static const int liquid_velocity_index = 0;
+    static const int pressure_index =
+        ShapeFunctionLiquidVelocity::NPOINTS * GlobalDim;
+
+    static const int liquid_velocity_size =
+        ShapeFunctionLiquidVelocity::NPOINTS * GlobalDim;
+    static const int pressure_size = ShapeFunctionPressure::NPOINTS;
+
+    using ShapeMatrixTypeLiquidVelocity =
+        ShapeMatrixPolicyType<ShapeFunctionLiquidVelocity, GlobalDim>;
+    using ShapeMatrixTypePressure =
+        ShapeMatrixPolicyType<ShapeFunctionPressure, GlobalDim>;
+
+    using NodalVectorType = typename ShapeMatrixTypePressure::NodalVectorType;
+
+public:
+    LocalAssemblerData(MeshLib::Element const& element,
+                       std::size_t const /*local_matrix_size*/,
+                       bool const is_axially_symmetric,
+                       unsigned const integration_order,
+                       StokesFlowProcessData const& process_data)
+        : _element(element),
+          _is_axially_symmetric(is_axially_symmetric),
+          _integration_method(integration_order),
+          _process_data(process_data)
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+        _ip_data.reserve(n_integration_points);
+
+        auto const shape_matrices_v =
+            NumLib::initShapeMatrices<ShapeFunctionLiquidVelocity,
+                                      ShapeMatrixTypeLiquidVelocity, GlobalDim>(
+                element, is_axially_symmetric, _integration_method);
+
+        auto const shape_matrices_p =
+            NumLib::initShapeMatrices<ShapeFunctionPressure,
+                                      ShapeMatrixTypePressure, GlobalDim>(
+                element, is_axially_symmetric, _integration_method);
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data.emplace_back(
+                shape_matrices_v[ip].N, shape_matrices_v[ip].dNdx,
+                shape_matrices_p[ip].N, shape_matrices_p[ip].dNdx,
+                _integration_method.getWeightedPoint(ip).getWeight() *
+                    shape_matrices_v[ip].integralMeasure *
+                    shape_matrices_v[ip].detJ);
+
+            auto& ip_data = _ip_data[ip];
+
+            ip_data.N_v_op = ShapeMatrixTypeLiquidVelocity::template MatrixType<
+                GlobalDim, liquid_velocity_size>::Zero(GlobalDim,
+                                                       liquid_velocity_size);
+
+            ip_data.dNdx_v_op =
+                ShapeMatrixTypeLiquidVelocity::template MatrixType<
+                    GlobalDim * GlobalDim,
+                    liquid_velocity_size>::Zero(GlobalDim * GlobalDim,
+                                                liquid_velocity_size);
+
+            for (int i = 0; i < GlobalDim; ++i)
+            {
+                ip_data.N_v_op
+                    .template block<1, liquid_velocity_size / GlobalDim>(
+                        i, i * liquid_velocity_size / GlobalDim)
+                    .noalias() = shape_matrices_v[ip].N;
+
+                ip_data.dNdx_v_op
+                    .template block<GlobalDim,
+                                    liquid_velocity_size / GlobalDim>(
+                        i * GlobalDim, i * liquid_velocity_size / GlobalDim)
+                    .noalias() = shape_matrices_v[ip].dNdx;
+            }
+        }
+    }
+
+    void assemble(double const t, double const dt,
+                  std::vector<double> const& local_x,
+                  std::vector<double> const& /*local_xdot*/,
+                  std::vector<double>& /*local_M_data*/,
+                  std::vector<double>& local_K_data,
+                  std::vector<double>& local_b_data) override
+    {
+        auto const local_matrix_size = liquid_velocity_size + pressure_size;
+        assert(local_x.size() == local_matrix_size);
+
+        auto local_p = Eigen::Map<const NodalVectorType>(
+            &local_x[pressure_index], pressure_size);
+
+        auto local_K = MathLib::createZeroedMatrix<
+            typename ShapeMatrixTypeLiquidVelocity::template MatrixType<
+                local_matrix_size, local_matrix_size>>(
+            local_K_data, local_matrix_size, local_matrix_size);
+        auto local_b = MathLib::createZeroedVector<
+            typename ShapeMatrixTypeLiquidVelocity::template VectorType<
+                local_matrix_size>>(local_b_data, local_matrix_size);
+
+        // Get block matrices
+        // Kvv
+        auto Kvv =
+            local_K.template block<liquid_velocity_size, liquid_velocity_size>(
+                liquid_velocity_index, liquid_velocity_index);
+        // Kvp
+        auto Kvp = local_K.template block<liquid_velocity_size, pressure_size>(
+            liquid_velocity_index, pressure_index);
+        // kpv
+        auto Kpv = local_K.template block<pressure_size, liquid_velocity_size>(
+            pressure_index, liquid_velocity_index);
+        // bv
+        auto bv = local_b.template segment<liquid_velocity_size>(
+            liquid_velocity_index);
+
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        ParameterLib::SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
+        auto const& b = _process_data.specific_body_force;
+
+        MaterialPropertyLib::VariableArray vars;
+
+        // create unit vector
+        assert(GlobalDim == 2);
+        auto const identity_mat =
+            Eigen::Matrix<double, GlobalDim, GlobalDim>::Identity().eval();
+        auto const unit_vector = Eigen::Map<const Eigen::RowVectorXd>(
+                                     identity_mat.data(), identity_mat.size())
+                                     .eval();
+
+        // Get material properties
+        auto const& medium =
+            *_process_data.media_map->getMedium(_element.getID());
+        auto const& phase = medium.phase("AqueousLiquid");
+
+        for (unsigned ip(0); ip < n_integration_points; ++ip)
+        {
+            pos.setIntegrationPoint(ip);
+
+            auto& ip_data = _ip_data[ip];
+
+            auto const& N_p = ip_data.N_p;
+            auto const& N_v_op = ip_data.N_v_op;
+
+            auto const& dNdx_v = ip_data.dNdx_v;
+            auto const& dNdx_v_op = ip_data.dNdx_v_op;
+            auto const& dNdx_p = ip_data.dNdx_p;
+
+            auto const& w = ip_data.integration_weight;
+
+            double p_int_pt = 0.0;
+
+            NumLib::shapeFunctionInterpolate(local_p, N_p, p_int_pt);
+
+            vars[static_cast<int>(
+                MaterialPropertyLib::Variable::phase_pressure)] = p_int_pt;
+
+            // Use the viscosity model to compute the viscosity
+            auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
+                                .template value<double>(vars, pos, t, dt);
+
+            // Kvv
+            for (int i = 0; i < GlobalDim; ++i)
+            {
+                Kvv.template block<liquid_velocity_size / GlobalDim,
+                                   liquid_velocity_size / GlobalDim>(
+                       i * liquid_velocity_size / GlobalDim,
+                       i * liquid_velocity_size / GlobalDim)
+                    .noalias() += w * dNdx_v.transpose() * mu * dNdx_v;
+            }
+
+            // Kvp
+            Kvp.noalias() += w * N_v_op.transpose() * dNdx_p;
+
+            // Kpv
+            Kpv.noalias() += w * N_p.transpose() * unit_vector * dNdx_v_op;
+
+            // bv
+            bv.noalias() += w * N_v_op.transpose() * b;
+        }
+    }
+
+    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+        const unsigned integration_point) const override
+    {
+        auto const& N_p = _ip_data[integration_point].N_p;
+
+        // assumes N_p is stored contiguously in memory
+        return Eigen::Map<const Eigen::RowVectorXd>(N_p.data(), N_p.size());
+    }
+
+    void computeSecondaryVariableConcrete(
+        double const /*t*/,
+        double const /*dt*/,
+        Eigen::VectorXd const& local_x,
+        Eigen::VectorXd const& /*local_x_dot*/) override
+    {
+        auto const local_p =
+            local_x.template segment<pressure_size>(pressure_index);
+
+        NumLib::interpolateToHigherOrderNodes<
+            ShapeFunctionPressure,
+            typename ShapeFunctionLiquidVelocity::MeshElement, GlobalDim>(
+            _element, _is_axially_symmetric, local_p,
+            *_process_data.pressure_interpolated);
+    }
+
+private:
+    MeshLib::Element const& _element;
+    bool const _is_axially_symmetric;
+    IntegrationMethod const _integration_method;
+    StokesFlowProcessData const& _process_data;
+
+    std::vector<IntegrationPointData<ShapeMatrixTypeLiquidVelocity,
+                                     ShapeMatrixTypePressure, GlobalDim,
+                                     ShapeFunctionLiquidVelocity::NPOINTS>,
+                Eigen::aligned_allocator<IntegrationPointData<
+                    ShapeMatrixTypeLiquidVelocity, ShapeMatrixTypePressure,
+                    GlobalDim, ShapeFunctionLiquidVelocity::NPOINTS>>>
+        _ip_data;
+};
+
+}  // namespace ProcessLib::StokesFlow