diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index cdbe4fd08a5c05cc652b1231c4911b8954219ffd..378af69166ee023088a2324acc9325dced3404f5 100644
--- a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
+++ b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
@@ -157,4 +157,16 @@ private:
     const MeshElementType* _ele;
 };
 
+/// Creates a TemplateIsoparametric element for the given shape functions and
+/// the underlying mesh element.
+template <typename ShapeFunction, typename ShapeMatricesType>
+NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>
+createIsoparametricFiniteElement(MeshLib::Element const& e)
+{
+    using FemType =
+        NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+
+    return FemType{
+        *static_cast<const typename ShapeFunction::MeshElement*>(&e)};
+}
 }  // namespace NumLib
diff --git a/NumLib/Function/Interpolation.h b/NumLib/Function/Interpolation.h
index f4a00c2f5ff3384f227346dc23d4a1ab66e274b9..64cc0d99d2c9af539a2cb27162d46cadb73b63a4 100644
--- a/NumLib/Function/Interpolation.h
+++ b/NumLib/Function/Interpolation.h
@@ -108,9 +108,10 @@ void interpolateToHigherOrderNodes(
     using SF = LowerOrderShapeFunction;
     using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
     using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
-    using FemType = TemplateIsoparametric<SF, ShapeMatricesType>;
 
-    FemType fe(*static_cast<const typename SF::MeshElement*>(&element));
+    auto const fe =
+        NumLib::createIsoparametricFiniteElement<SF, ShapeMatricesType>(
+            element);
     int const number_base_nodes = element.getNumberOfBaseNodes();
     int const number_all_nodes = element.getNumberOfNodes();
 
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
index 9c214e163e3e389e931bc09db99f37e8dd7330c4..254b0990350f7f0737898dff4aaff1667ba8ee9a 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -89,12 +89,8 @@ public:
     {
         (void)local_matrix_size; // unused, but needed for the interface
 
-
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
-        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &_surface_element));
+        auto const fe = NumLib::createIsoparametricFiniteElement<
+            ShapeFunction, ShapeMatricesType>(_surface_element);
 
         auto const n_integration_points =
             _integration_method.getNumberOfPoints();
diff --git a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
index 40437fa3317ff8863344972b7846068a8a1ccee6..c8d9dadc59efdfe03c2272ff9c2dc7822d8eb4a1 100644
--- a/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionLocalAssembler.h
@@ -52,10 +52,8 @@ public:
     {
         using ShapeMatricesType =
             ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &this->_element));
+        auto const fe = NumLib::createIsoparametricFiniteElement<
+            ShapeFunction, ShapeMatricesType>(Base::_element);
 
         unsigned const num_integration_points =
             Base::_integration_method.getNumberOfPoints();
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 3fd4bc8c49e4204ffed3e8af072012696d54acc7..26320aaa0aaf4862cdf19e19ada7dec3ffc864ea 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -809,11 +809,8 @@ public:
     {
         // eval shape matrices at given point
         auto const shape_matrices = [&]() {
-            using FemType =
-                NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
-            FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-                &_element));
+            auto const fe = NumLib::createIsoparametricFiniteElement<
+                ShapeFunction, ShapeMatricesType>(_element);
 
             typename ShapeMatricesType::ShapeMatrices shape_matrices(
                 ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index c4be98e8eacbe3ef909294170144b6c07e9b0149..6d14083c922ff47e7bcf334fed2cd1d24262b023 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -138,11 +138,8 @@ public:
                             std::vector<double> const& local_x) const override
     {
         // eval dNdx and invJ at p
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
-        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &_element));
+        auto const fe = NumLib::createIsoparametricFiniteElement<
+            ShapeFunction, ShapeMatricesType>(_element);
 
         typename ShapeMatricesType::ShapeMatrices shape_matrices(
             ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index 55bc5c3957574713bffb18ca4919e853b092d94c..6475ea65d75bbd55ede2e40f5a85f9178558076b 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -132,11 +132,8 @@ public:
                             std::vector<double> const& local_x) const override
     {
         // eval dNdx and invJ at given point
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
-        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &this->_element));
+        auto const fe = NumLib::createIsoparametricFiniteElement<
+            ShapeFunction, ShapeMatricesType>(_element);
 
         typename ShapeMatricesType::ShapeMatrices shape_matrices(
             ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
index 3e609c48ef775f1f6ef64cdfca2b33ef414f1cea..99b41376e121c4a8b735d12cde4e90a9b097ed01 100644
--- a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
@@ -92,10 +92,8 @@ public:
     {
         using ShapeMatricesType =
             ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &_element));
+        auto const fe = NumLib::createIsoparametricFiniteElement<
+            ShapeFunction, ShapeMatricesType>(_element);
 
         unsigned const num_integration_points =
             _integration_method.getNumberOfPoints();
diff --git a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
index e281e7be17999653ad1e798ec8aa8b273c78f754..0829219f32d1254d1c02031f5a588dfe5d59c879 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
@@ -71,11 +71,8 @@ public:
           _bulk_element_id(bulk_element_ids[surface_element.getID()]),
           _bulk_face_id(bulk_face_ids[surface_element.getID()])
     {
-        using FemType =
-            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
-        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
-            &surface_element));
+        auto const fe = NumLib::createIsoparametricFiniteElement<
+            ShapeFunction, ShapeMatricesType>(_surface_element);
 
         std::size_t const n_integration_points =
             _integration_method.getNumberOfPoints();
diff --git a/ProcessLib/Utils/InitShapeMatrices.h b/ProcessLib/Utils/InitShapeMatrices.h
index 773188972bbb5770a8e984a0ae984f5da6f3e5ee..a038837a7380ee4fd33e06b800d3303782f678cc 100644
--- a/ProcessLib/Utils/InitShapeMatrices.h
+++ b/ProcessLib/Utils/InitShapeMatrices.h
@@ -28,10 +28,9 @@ initShapeMatrices(MeshLib::Element const& e, bool is_axially_symmetric,
         Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
         shape_matrices;
 
-    using FemType = NumLib::TemplateIsoparametric<
-        ShapeFunction, ShapeMatricesType>;
-
-    FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(&e));
+    auto const fe =
+        NumLib::createIsoparametricFiniteElement<ShapeFunction,
+                                                 ShapeMatricesType>(e);
 
     unsigned const n_integration_points = integration_method.getNumberOfPoints();
 
@@ -52,10 +51,9 @@ double interpolateXCoordinate(
     MeshLib::Element const& e,
     typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
 {
-    using FemType = NumLib::TemplateIsoparametric<
-        ShapeFunction, ShapeMatricesType>;
-
-    FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(&e));
+    auto const fe =
+        NumLib::createIsoparametricFiniteElement<ShapeFunction,
+                                                 ShapeMatricesType>(e);
 
     return fe.interpolateZerothCoordinate(N);
 }
@@ -65,10 +63,9 @@ std::array<double, 3> interpolateCoordinates(
     MeshLib::Element const& e,
     typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
 {
-    using FemType = NumLib::TemplateIsoparametric<
-        ShapeFunction, ShapeMatricesType>;
-
-    FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(&e));
+    auto const fe =
+        NumLib::createIsoparametricFiniteElement<ShapeFunction,
+                                                 ShapeMatricesType>(e);
 
     return fe.interpolateCoordinates(N);
 }
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
index ae5181077468a2cf14197397bf0050156c05eb28..43004901aac26f943ab9edcb7a8d4f11fdcb9147 100644
--- a/Tests/NumLib/TestFunctionInterpolation.cpp
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -47,8 +47,6 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
     using ShapeFunction = NumLib::ShapeLine2;
     using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, 1u>;
 
-    using FemType = NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
-
     using IntegrationMethod = NumLib::GaussLegendreIntegrationPolicy<
         ShapeFunction::MeshElement>::IntegrationMethod;
 
@@ -60,7 +58,9 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
         MeshLib::Line(std::array<MeshLib::Node*, 2>{{&pt_a, &pt_b}});
 
     // set up shape function
-    FemType finite_element(*static_cast<const ShapeFunction::MeshElement*>(&element));
+    auto const finite_element =
+        NumLib::createIsoparametricFiniteElement<ShapeFunction,
+                                                 ShapeMatricesType>(element);
 
     const unsigned integration_order = 2;
     IntegrationMethod integration_method(integration_order);