Skip to content

Commit

Permalink
Replaced initialization and obvious computation with STL algorithms
Browse files Browse the repository at this point in the history
  • Loading branch information
JohanMabille committed Sep 25, 2020
1 parent 6307c44 commit 9a49a6c
Showing 1 changed file with 44 additions and 70 deletions.
114 changes: 44 additions & 70 deletions proteus/mprans/RANS2P2D.h
Original file line number Diff line number Diff line change
Expand Up @@ -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>;

Expand Down Expand Up @@ -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++)
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down

0 comments on commit 9a49a6c

Please sign in to comment.