From b63c3840cc2f1dd7d20cb0f3306489dfc88015bf Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Mon, 20 May 2019 14:58:01 +0200
Subject: [PATCH] [NL] Add a create function for common FE type.

---
 .../Fem/FiniteElement/TemplateIsoparametric.h | 12 +++++++++++
 NumLib/Function/Interpolation.h               |  5 +++--
 ...DirichletBoundaryConditionLocalAssembler.h |  8 ++-----
 .../PythonBoundaryConditionLocalAssembler.h   |  6 ++----
 .../ComponentTransportFEM.h                   |  7 ++-----
 .../GroundwaterFlow/GroundwaterFlowFEM.h      |  7 ++-----
 ProcessLib/HT/HTFEM.h                         |  7 ++-----
 .../Python/PythonSourceTermLocalAssembler.h   |  6 ++----
 .../SurfaceFlux/SurfaceFluxLocalAssembler.h   |  7 ++-----
 ProcessLib/Utils/InitShapeMatrices.h          | 21 ++++++++-----------
 Tests/NumLib/TestFunctionInterpolation.cpp    |  6 +++---
 11 files changed, 41 insertions(+), 51 deletions(-)

diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index cdbe4fd08a5..378af69166e 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 f4a00c2f5ff..64cc0d99d2c 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 9c214e163e3..254b0990350 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 40437fa3317..c8d9dadc59e 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 3fd4bc8c49e..26320aaa0aa 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 c4be98e8eac..6d14083c922 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 55bc5c39575..6475ea65d75 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 3e609c48ef7..99b41376e12 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 e281e7be179..0829219f32d 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 773188972bb..a038837a738 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 ae518107746..43004901aac 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);
-- 
GitLab