[PL/HT] Use shape matrix cache in HT assembler

ProcessLib/HT/HTFEM.h

@@ -91,7 +91,8 @@ public:
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override
     {
-        auto const& N = _ip_data[integration_point].N;
+        auto const& N = _shape_matrix_cache.NsHigherOrder<
+            typename ShapeFunction::MeshElement>()[integration_point];
 
         // assumes N is stored contiguously in memory
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
@@ -273,11 +274,15 @@ protected:
             *_process_data.media_map.getMedium(_element.getID());
         auto const& liquid_phase = medium.phase("AqueousLiquid");
 
+        auto const& Ns =
+            _shape_matrix_cache
+                .NsHigherOrder<typename ShapeFunction::MeshElement>();
+
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
             auto const& ip_data = _ip_data[ip];
-            auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
+            auto const& N = Ns[ip];
 
             pos.setIntegrationPoint(ip);
 

ProcessLib/HT/MonolithicHTFEM.h

@@ -123,13 +123,17 @@ public:
         double average_velocity_norm = 0.0;
         ip_flux_vector.reserve(n_integration_points);
 
+        auto const& Ns =
+            this->_shape_matrix_cache
+                .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
         for (unsigned ip(0); ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
 
             auto const& ip_data = this->_ip_data[ip];
-            auto const& N = ip_data.N;
             auto const& dNdx = ip_data.dNdx;
+            auto const& N = Ns[ip];
             auto const& w = ip_data.integration_weight;
 
             double T_int_pt = 0.0;

ProcessLib/HT/StaggeredHTFEM-impl.h

@@ -79,13 +79,17 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHydraulicEquation(
     unsigned const n_integration_points =
         this->_integration_method.getNumberOfPoints();
 
+    auto const& Ns =
+        this->_shape_matrix_cache
+            .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
     for (unsigned ip(0); ip < n_integration_points; ip++)
     {
         pos.setIntegrationPoint(ip);
 
         auto const& ip_data = this->_ip_data[ip];
-        auto const& N = ip_data.N;
         auto const& dNdx = ip_data.dNdx;
+        auto const& N = Ns[ip];
         auto const& w = ip_data.integration_weight;
 
         double p_int_pt = 0.0;
@@ -208,13 +212,17 @@ void StaggeredHTFEM<ShapeFunction, GlobalDim>::assembleHeatTransportEquation(
     double average_velocity_norm = 0.0;
     ip_flux_vector.reserve(n_integration_points);
 
+    auto const& Ns =
+        this->_shape_matrix_cache
+            .template NsHigherOrder<typename ShapeFunction::MeshElement>();
+
     for (unsigned ip(0); ip < n_integration_points; ip++)
     {
         pos.setIntegrationPoint(ip);
 
         auto const& ip_data = this->_ip_data[ip];
-        auto const& N = ip_data.N;
         auto const& dNdx = ip_data.dNdx;
+        auto const& N = Ns[ip];
         auto const& w = ip_data.integration_weight;
 
         double p_at_xi = 0.;