From 9a49a6c59955850343e1a0f7c7d8883ededf78f7 Mon Sep 17 00:00:00 2001
From: Johan Mabille <johan.mabille@gmail.com>
Date: Thu, 24 Sep 2020 17:09:59 +0200
Subject: [PATCH] Replaced initialization and obvious computation with STL
 algorithms

---
 proteus/mprans/RANS2P2D.h | 114 +++++++++++++++-----------------------
 1 file changed, 44 insertions(+), 70 deletions(-)

diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h
index a490ad70df..6bc5711ce3 100644
--- a/proteus/mprans/RANS2P2D.h
+++ b/proteus/mprans/RANS2P2D.h
@@ -24,6 +24,17 @@ const  double DM3=1.0;//1-point-wise divergence, 0-point-wise rate of volume cha
 const double inertial_term=1.0;
 namespace proteus
 {
+    
+    namespace detail
+    {
+        template <size_t N, class T>
+        inline std::array<T, N> make_array(const T& t)
+        {
+            std::array<T, N> res;
+            std::fill_n(res.begin(), N, t);
+        }
+    }
+
   template<int nSpace, int nP, int nQ, int nEBQ>
   using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nQ, nEBQ>;
 
@@ -1682,35 +1693,22 @@ namespace proteus
       for(int eN=0;eN<nElements_global;eN++)
         {
           //declare local storage for element residual and initialize
-          register double elementResidual_p[nDOF_test_element],elementResidual_p_check[nDOF_test_element],elementResidual_mesh[nDOF_test_element],
-            elementResidual_u[nDOF_v_test_element],
-            elementResidual_v[nDOF_v_test_element],
-            pelementResidual_u[nDOF_v_test_element],
-            pelementResidual_v[nDOF_v_test_element],
-            velocityErrorElement[nDOF_v_test_element],
-            eps_rho,eps_mu;
+
+          register auto elementResidual_p       = detail::make_array<nDOF_test_element>(0.0);
+          register auto elementResidual_p_check = detail::make_array<nDOF_test_element>(0.0);
+          register auto elementResidual_mesh    = detail::make_array<nDOF_test_element>(0.0);
+          register auto elementResidual_u       = detail::make_array<nDOF_v_test_element>(0.0);
+          register auto elementResidual_v       = detail::make_array<nDOF_v_test_element>(0.0);
+          register auto pelementResidual_u      = detail::make_array<nDOF_v_test_element>(0.0);
+          register auto pelementResidual_v      = detail::make_array<nDOF_v_test_element>(0.0);
+          register auto velocityErrorElement    = detail::make_array<nDOF_v_test_element>(0.0);
+          double eps_rho,eps_mu;
           bool element_active=false;
           elementIsActive[eN]=false;
           const double* elementResidual_w(NULL);
           double mesh_volume_conservation_element=0.0,
             mesh_volume_conservation_element_weak=0.0;
-          for (int i=0;i<nDOF_test_element;i++)
-            {
-              int eN_i = eN*nDOF_test_element+i;
-              elementResidual_p_save.data()[eN_i]=0.0;
-              elementResidual_mesh[i]=0.0;
-              elementResidual_p[i]=0.0;
-              elementResidual_p_check[i]=0.0;
-            }
-          for (int i=0;i<nDOF_v_test_element;i++)
-            {
-              elementResidual_u[i]=0.0;
-              elementResidual_v[i]=0.0;
-              pelementResidual_u[i]=0.0;
-              pelementResidual_v[i]=0.0;
-              velocityErrorElement[i]=0.0;
-            }//i
-          //Use for plotting result
+          std::fill_n(elementResidual_p_save.data() + eN*nDOF_test_element, nDOF_test_element, 0.0);
           if(use_ball_as_particle==1 && nParticles > 0)
             {
               for (int I=0;I<nDOF_mesh_trial_element;I++)
@@ -1925,6 +1923,7 @@ namespace proteus
                     dmass_ham_u_s=0.0,
                     dmass_ham_v_s=0.0,
                     dmass_ham_w_s=0.0;
+
                   //get jacobian, etc for mapping reference element
                   gf_s.set_quad(k);
                   gf.set_quad(k);
@@ -1965,18 +1964,12 @@ namespace proteus
                   //get the trial function gradients
                   ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
                   ck_v.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_v_trial_element*nSpace],jacInv,vel_grad_trial);
-                  for (int i=0; i < nDOF_trial_element; i++)
-                    {
-                      p_trial[i] = p_trial_ref.data()[k*nDOF_trial_element + i];
-                      p_grad_trial_ib[i*nSpace + 0] = p_grad_trial[i*nSpace+0];
-                      p_grad_trial_ib[i*nSpace + 1] = p_grad_trial[i*nSpace+1];
-                    }
-                  for (int i=0; i < nDOF_v_trial_element; i++)
-                    {
-                      vel_trial[i] = vel_trial_ref.data()[k*nDOF_v_trial_element + i];
-                      vel_grad_trial_ib[i*nSpace + 0] = vel_grad_trial[i*nSpace+0];
-                      vel_grad_trial_ib[i*nSpace + 1] = vel_grad_trial[i*nSpace+1];
-                    }
+
+                  std::copy_n(p_trial_ref.data() + k*nDOF_trial_element, nDOF_trial_element, p_trial);
+                  std::copy(std::begin(p_grad_trial), std::end(p_grad_trial), std::begin(p_grad_trial_ib));
+                  std::copy_n(vel_trial_ref.data() + k*nDOF_v_trial_element, nDOF_v_trial_element, p_trial);
+                  std::copy(std::begin(vel_grad_trial), std::end(vel_grad_trial), std::begin(vel_grad_trial_ib));
+
                   if (icase == 0)
                     {
 #ifdef IFEMBASIS
@@ -2092,42 +2085,23 @@ namespace proteus
                   if (PRESSURE_PROJECTION_STABILIZATION)
                     ck.DOFaverage(p_dof.data(), &p_l2g.data()[eN_nDOF_trial_element],p_element_avg);
                   //precalculate test function products with integration weights
+                  auto lambda_dv = [&dV](double d) { return d * dV; };
 #ifdef IFEMGALERKIN
-                  for (int j=0;j<nDOF_test_element;j++)
-                    {
-                      p_test_dV[j] = p_trial[j]*dV;
-                      for (int I=0;I<nSpace;I++)
-                        {
-                          p_grad_test_dV[j*nSpace+I]   = p_grad_trial_ib[j*nSpace+I]*dV;
-                        }
-                    }
-                  //precalculate test function products with integration weights
-                  for (int j=0;j<nDOF_v_test_element;j++)
-                    {
-                      vel_test_dV[j] = vel_trial[j]*dV;
-                      for (int I=0;I<nSpace;I++)
-                        {
-                          vel_grad_test_dV[j*nSpace+I] = vel_grad_trial_ib[j*nSpace+I]*dV;
-                        }
-                    }
+                  std::transform(std::begin(p_trial), std::end(p_trial), std::begin(p_test_dV), lambda_dv);
+                  std::transform(std::begin(p_grad_trial_ib), std::end(p_grad_trial_ib),
+                                 std::begin(p_grad_test_dV), lambda_dv);
+                  std::transform(std::begin(vel_trial), std::end(vel_trial), std::begin(vel_test_dV), lambda_dv);
+                  std::transform(std::begin(vel_grad_trial_ib), std::end(vel_grad_trial_ib),
+                                 std::begin(vel_grad_test_dV), lambda_dv);
 #else
-                  for (int j=0;j<nDOF_test_element;j++)
-                    {
-                      p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
-                      for (int I=0;I<nSpace;I++)
-                        {
-                          p_grad_test_dV[j*nSpace+I]   = p_grad_trial[j*nSpace+I]*dV;//assume test_i = trial_i, not using ib basis here
-                        }
-                    }
-                  //precalculate test function products with integration weights
-                  for (int j=0;j<nDOF_v_test_element;j++)
-                    {
-                      vel_test_dV[j] = vel_test_ref.data()[k*nDOF_v_trial_element+j]*dV;
-                      for (int I=0;I<nSpace;I++)
-                        {
-                          vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;//assume test_i = trial_i
-                        }
-                    }
+                  double* p_test_ref_beg = p_test_ref.data() + k*nDOF_trial_element;
+                  std::transform(p_test_ref_beg, p_test_ref_beg+nDOF_test_element, p_test_dV, lambda_dv);
+                  std::transform(std::begin(p_grad_trial), std::end(p_grad_trial),
+                                 std::begin(p_grad_test_dV), lambda_dv);
+                  double* vel_test_ref_beg = vel_test_ref.data() + k*nDOF_v_trial_element;
+                  std::transform(vel_test_ref_beg, vel_test_ref_beg+nDOF_v_test_element, vel_test_dV, lambda_dv);
+                  std::transform(std::begin(vel_grad_trial), std::end(vel_grad_trial),
+                                 std::begin(vel_grad_test_dV), lambda_dv);
 #endif
                   //todo: extend this to higher-order meshes, for now assume mesh trial and p trial are same 
                   double div_mesh_velocity=0.0;