diff --git a/ProcessLib/StokesFlow/CreateLocalAssemblers.h b/ProcessLib/StokesFlow/CreateLocalAssemblers.h
new file mode 100644
index 0000000000000000000000000000000000000000..9d3215354d370ba2a6e5761f299decd0d0b7c0fe
--- /dev/null
+++ b/ProcessLib/StokesFlow/CreateLocalAssemblers.h
@@ -0,0 +1,87 @@
+/**
+ * \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
diff --git a/ProcessLib/StokesFlow/LocalDataInitializer.h b/ProcessLib/StokesFlow/LocalDataInitializer.h
new file mode 100644
index 0000000000000000000000000000000000000000..caaaee96fce7b5f40cacd9a8ff2ca76c1612aa53
--- /dev/null
+++ b/ProcessLib/StokesFlow/LocalDataInitializer.h
@@ -0,0 +1,287 @@
+/**
+ * \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
diff --git a/ProcessLib/StokesFlow/StokesFlowFEM.h b/ProcessLib/StokesFlow/StokesFlowFEM.h
new file mode 100644
index 0000000000000000000000000000000000000000..cd3a9f3ab4611e739f494dca071faea957230c5f
--- /dev/null
+++ b/ProcessLib/StokesFlow/StokesFlowFEM.h
@@ -0,0 +1,257 @@
+/**
+ * \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