From 5ed03935142e96085ee5d00b442010e3f3d8644b Mon Sep 17 00:00:00 2001
From: Johan Mabille <johan.mabille@gmail.com>
Date: Thu, 24 Sep 2020 11:03:32 +0200
Subject: [PATCH] Replaced initialization loops and trivial computations with
 STL algorithms

---
 proteus/mprans/ArgumentsDict.cpp |   2 +-
 proteus/mprans/RANS2P2D.h        | 114 ++++++++++++-------------------
 2 files changed, 43 insertions(+), 73 deletions(-)

diff --git a/proteus/mprans/ArgumentsDict.cpp b/proteus/mprans/ArgumentsDict.cpp
index 0934af1b98..655e2c2a3e 100644
--- a/proteus/mprans/ArgumentsDict.cpp
+++ b/proteus/mprans/ArgumentsDict.cpp
@@ -119,7 +119,7 @@ namespace proteus
                 });
         // IMPORTANT: keep the scalar overloads AFTER the pyarray overloads,
         // because numpy array containing a single elements can be implicitly
-        // converted to scalars. If the scalar overload is define before the
+        // converted to scalars. If the scalar overload is defined before the
         // pyarray overload, it is chosen by pybind11 and the array ends up
         // in the wrong dictionary.
         cl.def("__setitem__",
diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h
index a490ad70df..af2a8af139 100644
--- a/proteus/mprans/RANS2P2D.h
+++ b/proteus/mprans/RANS2P2D.h
@@ -24,6 +24,16 @@ 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>;
 
@@ -1617,7 +1627,7 @@ namespace proteus
       xt::pyarray<double>& forcex = args.array<double>("forcex");
       xt::pyarray<double>& forcey = args.array<double>("forcey");
       xt::pyarray<double>& forcez = args.array<double>("forcez");
-      int     use_ball_as_particle = args.scalar<int>("use_ball_as_particle");
+      int use_ball_as_particle = args.scalar<int>("use_ball_as_particle");
       xt::pyarray<double>& ball_center = args.array<double>("ball_center");
       xt::pyarray<double>& ball_radius = args.array<double>("ball_radius");
       xt::pyarray<double>& ball_velocity = args.array<double>("ball_velocity");
@@ -1680,36 +1690,21 @@ namespace proteus
       /*          <<"Ball Info: velocity "<<ball_velocity[0]<<'\t'<<ball_velocity[1]<<'\t'<<ball_velocity[2]<<std::endl */
       /*          <<"Ball Info: angular "<<ball_angular_velocity[0]<<ball_angular_velocity[1]<<ball_angular_velocity[2]<<std::endl; */
       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;
+      {
+          auto elementResidual_p       = detail::make_array<nDOF_test_element>(0.0);
+          auto elementResidual_p_check = detail::make_array<nDOF_test_element>(0.0);
+          auto elementResidual_mesh    = detail::make_array<nDOF_test_element>(0.0);
+          auto elementResidual_u       = detail::make_array<nDOF_v_test_element>(0.0);
+          auto elementResidual_v       = detail::make_array<nDOF_v_test_element>(0.0);
+          auto pelementResidual_u      = detail::make_array<nDOF_v_test_element>(0.0);
+          auto pelementResidual_v      = detail::make_array<nDOF_v_test_element>(0.0);
+          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
           if(use_ball_as_particle==1 && nParticles > 0)
             {
@@ -1925,6 +1920,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 +1961,11 @@ 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_n(p_grad_trial, nDOF_trial_element*nSpace, p_grad_trial_ib);
+                  std::copy_n(vel_trial_ref.data() + k*nDOF_v_trial_element, nDOF_v_trial_element, vel_trial);
+                  std::copy_n(vel_grad_trial, nDOF_v_trial_element*nSpace, vel_grad_trial_ib);
+
                   if (icase == 0)
                     {
 #ifdef IFEMBASIS
@@ -2092,42 +2081,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 dv_lambda = [dV](double arg) { return arg * 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(p_trial, p_trial + nDOF_test_element, p_test_dV, dv_lambda)
+                  std::transform(p_grad_trial_ib, p_grad_trial_ib + nDOF_test_element*nSpace, p_grad_test_dV, dv_lambda);
+                  std::transform(vel_trial, vel_trial + nDOF_v_test_element, vel_test_dV, dv_lambda);
+                  std::transform(vel_grad_trial_ib, vel_grad_trial_ib + nDOF_v_test_element*nSpace, vel_grad_test_dV, dv_lambda);
 #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
-                        }
-                    }
+                  std::transform(p_test_ref.data() + k*nDOF_trial_element,
+                                 p_test_ref.data() + (k+1)*nDOF_trial_element,
+                                 p_test_dV,
+                                 dv_lambda);
+                  std::transform(p_grad_trial, p_grad_trial+nDOF_test_element*nSpace, p_grad_test_dV, dv_lambda);
+                  std::transform(vel_test_ref.data() + k*nDOF_v_trial_element,
+                                 vel_test_ref.data() + (k+1)*nDOF_v_test_element,
+                                 vel_test_dV,
+                                 dv_lambda);
+                  std::transform(vel_grad_trial, vel_grad_trial+nDOF_v_test_element*nSpace, vel_grad_test_dV, dv_lambda);
 #endif
                   //todo: extend this to higher-order meshes, for now assume mesh trial and p trial are same 
                   double div_mesh_velocity=0.0;