Skip to content

Commit

Permalink
Replaced initialization loops and trivial computations with STL algor…
Browse files Browse the repository at this point in the history 
…ithms
  • Loading branch information
@JohanMabille
JohanMabille committed Sep 24, 2020
1 parent 6307c44 commit 5ed0393
Show file tree
Hide file tree
Showing 2 changed files with 43 additions and 73 deletions.
2 changes: 1 addition & 1 deletion proteus/mprans/ArgumentsDict.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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__",
Expand Down
114 changes: 42 additions & 72 deletions proteus/mprans/RANS2P2D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>;

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

0 comments on commit 5ed0393

Please sign in to comment.