diff --git a/proteus/mprans/RANS2P2D.h b/proteus/mprans/RANS2P2D.h index da130b6cf6..c4819ba865 100644 --- a/proteus/mprans/RANS2P2D.h +++ b/proteus/mprans/RANS2P2D.h @@ -17,6 +17,8 @@ #include "mpi.h" #include "proteus_lapack.h" +#include + namespace py = pybind11; #define ZEROVEC {0.,0.} @@ -28,6 +30,19 @@ 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 + inline std::array make_array(const T& t) + { + std::array res; + std::fill_n(res.begin(), N, t); + return res; + } + } + inline double enorm(double* v) { return std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); @@ -118,6 +133,7 @@ namespace proteus } } + template using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix; @@ -1942,18 +1958,24 @@ namespace proteus for(int eN=0;eN(0.0); + auto elementResidual_p_check = detail::make_array(0.0); + auto elementResidual_mesh = detail::make_array(0.0); + auto elementResidual_u = detail::make_array(0.0); + auto elementResidual_v = detail::make_array(0.0); + auto pelementResidual_u = detail::make_array(0.0); + auto pelementResidual_v = detail::make_array(0.0); + auto velocityErrorElement = detail::make_array(0.0); + double eps_rho,eps_mu; bool element_active=false; isActiveElement[eN]=0; const double* elementResidual_w(NULL); double mesh_volume_conservation_element=0.0, mesh_volume_conservation_element_weak=0.0; + + std::fill_n(elementResidual_p_save.data() + eN*nDOF_test_element, nDOF_test_element, 0.0); + int particle_index=0; for (int i=0;i 0) { double min_d = 1e10; for (int I=0;I 1.0e-8) + if (fabs(p_trial_ref(k, vi) - p_trial[vi]) > 1.0e-8) { for (int vj=0; vj < nDOF_trial_element; vj++) - std::cout<<"Trial "< 1.0e-8) @@ -2329,10 +2345,10 @@ namespace proteus prob=true; } //velocity - if (fabs(vel_trial_ref.data()[k*nDOF_v_trial_element + vi] - vel_trial[vi]) > 1.0e-8) + if (fabs(vel_trial_ref(k, vi) - vel_trial[vi]) > 1.0e-8) { for (int vj=0; vj < nDOF_v_trial_element; vj++) - std::cout<<"Trial "< 1.0e-8) @@ -2392,7 +2408,7 @@ namespace proteus #else for (int j=0;j 0) { - int ball_index=get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x,y,z,distance_to_solids.data()[eN_k]); + int ball_index=get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x,y,z,distance_to_solids(eN, k)); get_normal_to_ith_ball(nParticles, ball_center.data(), ball_radius.data(),ball_index,x,y,z,ball_n[0],ball_n[1]); } else @@ -2437,25 +2453,19 @@ namespace proteus //distance_to_solids is given in Prestep } if (nParticles > 0) - phi_solid.data()[eN_k] = distance_to_solids.data()[eN_k]; - const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]); + phi_solid(eN, k) = distance_to_solids(eN, k); + const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter(eN)); // //calculate pde coefficients at quadrature points // - const double H_s = gf_s.H(particle_eps,phi_solid.data()[eN_k]); - const double D_s = gf_s.D(particle_eps,phi_solid.data()[eN_k]); + const double H_s = gf_s.H(particle_eps,phi_solid(eN, k)); + const double D_s = gf_s.D(particle_eps,phi_solid(eN, k)); //save velocity at quadrature points for other models to use - double p_e = q_u_0.data()[eN_k] - p, - u_e = q_u_1.data()[eN_k] - u, - v_e = q_u_2.data()[eN_k] - v, + double p_e = q_u_0(eN ,k) - p, + u_e = q_u_1(eN ,k) - u, + v_e = q_u_2(eN, k) - v, velocity_e=sqrt(u_e*u_e + v_e*v_e); - - /* q_u_0.data()[eN_k] = p; */ - /* q_u_1.data()[eN_k] = u; */ - /* q_u_2.data()[eN_k] = v; */ - /* q_u_3.data()[eN_k] = 0.0; */ - double rho,nu; if (gf.useExact) { @@ -2497,18 +2507,18 @@ namespace proteus sigma, rho, nu, - elementDiameter.data()[eN], + elementDiameter(eN), smagorinskyConstant, turbulenceClosureModel, g.data(), useVF, - vf.data()[eN_k], - phi.data()[eN_k], - &normal_phi.data()[eN_k_nSpace], - kappa_phi.data()[eN_k], + vf(eN, k), + phi(eN, k), + &normal_phi(eN, k, 0), + kappa_phi(eN, k), //VRANS porosity, - phi_solid.data()[eN_k],//distance to solid + phi_solid(eN, k),//distance to solid p_old, u_old, v_old, @@ -2527,8 +2537,8 @@ namespace proteus v, w, LAG_LES, - q_eddy_viscosity.data()[eN_k], - q_eddy_viscosity_last.data()[eN_k], + q_eddy_viscosity(eN, k), + q_eddy_viscosity_last(eN, k), mom_u_acc, dmom_u_acc_u, mom_v_acc, @@ -2581,20 +2591,20 @@ namespace proteus dmom_w_ham_u, dmom_w_ham_v, dmom_w_ham_w, - forcex.data()[eN_k], - forcey.data()[eN_k], - forcez.data()[eN_k]); - q_rho.data()[eN_k] = rho; + forcex(eN, k), + forcey(eN, k), + forcez(eN, k)); + q_rho(eN, k) = rho; //VRANS - mass_source = q_mass_source.data()[eN_k]; + mass_source = q_mass_source(eN, k); //todo: decide if these should be lagged or not? updateDarcyForchheimerTerms_Ergun(NONCONSERVATIVE_FORM, /* linearDragFactor, */ /* nonlinearDragFactor, */ /* porosity, */ /* meanGrainSize, */ - q_dragAlpha.data()[eN_k], - q_dragBeta.data()[eN_k], + q_dragAlpha(eN, k), + q_dragBeta(eN, k), eps_rho, eps_mu, rho_0, @@ -2602,19 +2612,19 @@ namespace proteus rho_1, nu_1, useVF, - vf.data()[eN_k], - phi.data()[eN_k], + vf(eN, k), + phi(eN, k), u, v, w, - q_velocity_sge.data()[eN_k_nSpace+0], - q_velocity_sge.data()[eN_k_nSpace+1], - q_velocity_sge.data()[eN_k_nSpace+1],//dummy entry for 2D - eps_porous.data()[elementFlags.data()[eN]], - phi_porous.data()[eN_k], - q_velocity_porous.data()[eN_k_nSpace+0], - q_velocity_porous.data()[eN_k_nSpace+1], - q_velocity_porous.data()[eN_k_nSpace+1],//dummy entry for 2D + q_velocity_sge(eN, k, 0), + q_velocity_sge(eN, k, 1), + q_velocity_sge(eN, k, 1),//dummy entry for 2D + eps_porous(elementFlags(eN)), + phi_porous(eN, k), + q_velocity_porous(eN, k, 0), + q_velocity_porous(eN, k, 1), + q_velocity_porous(eN, k, 1),//dummy entry for 2D mom_u_source, mom_v_source, mom_w_source, @@ -2634,14 +2644,14 @@ namespace proteus rho_1, nu_1, useVF, - vf.data()[eN_k], - phi.data()[eN_k], + vf(eN, k), + phi(eN, k), porosity, c_mu, //mwf hack - q_turb_var_0.data()[eN_k], - q_turb_var_1.data()[eN_k], - &q_turb_var_grad_0.data()[eN_k_nSpace], - q_eddy_viscosity.data()[eN_k], + q_turb_var_0(eN, k), + q_turb_var_1(eN, k), + &q_turb_var_grad_0(eN, k, 0), + q_eddy_viscosity(eN, k), mom_uu_diff_ten, mom_vv_diff_ten, mom_ww_diff_ten, @@ -2689,17 +2699,17 @@ namespace proteus // //calculate time derivative at quadrature points // - if (q_dV_last.data()[eN_k] <= -100) - q_dV_last.data()[eN_k] = dV; - q_dV.data()[eN_k] = dV; + if (q_dV_last(eN, k) <= -100) + q_dV_last(eN, k) = dV; + q_dV(eN, k) = dV; ck.bdf(alphaBDF, - q_mom_u_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV, + q_mom_u_acc_beta_bdf(eN, k)*q_dV_last(eN, k)/dV, mom_u_acc, dmom_u_acc_u, mom_u_acc_t, dmom_u_acc_u_t); ck.bdf(alphaBDF, - q_mom_v_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV, + q_mom_v_acc_beta_bdf(eN, k)*q_dV_last(eN, k)/dV, mom_v_acc, dmom_v_acc_v, mom_v_acc_t, @@ -2716,20 +2726,20 @@ namespace proteus //calculate strong residual pdeResidual_p = ck.Advection_strong(dmass_adv_u,grad_u) + ck.Advection_strong(dmass_adv_v,grad_v) + - DM2*MOVING_DOMAIN*ck.Reaction_strong(alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity) + + DM2*MOVING_DOMAIN*ck.Reaction_strong(alphaBDF*(dV-q_dV_last(eN, k))/dV - div_mesh_velocity) + ck.Reaction_strong(mass_source); if (NONCONSERVATIVE_FORM > 0.0) { dmom_adv_sge[0] = 0.0; dmom_adv_sge[1] = 0.0; - dmom_ham_grad_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt); - dmom_ham_grad_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt); + dmom_ham_grad_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 0) - MOVING_DOMAIN*xt); + dmom_ham_grad_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 1) - MOVING_DOMAIN*yt); } else { - dmom_adv_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt); - dmom_adv_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt); + dmom_adv_sge[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 0) - MOVING_DOMAIN*xt); + dmom_adv_sge[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 1) - MOVING_DOMAIN*yt); dmom_ham_grad_sge[0] = 0.0; dmom_ham_grad_sge[1] = 0.0; } @@ -2755,7 +2765,7 @@ namespace proteus //add contributions from mass and source terms double tmpR=dmom_u_acc_u_t + dmom_u_source[0]; calculateSubgridError_tau(hFactor, - elementDiameter.data()[eN], + elementDiameter(eN), tmpR,//dmom_u_acc_u_t, dmom_u_acc_u, mv_tau,//dmom_adv_sge, @@ -2763,7 +2773,7 @@ namespace proteus dmom_u_ham_grad_p[0], tau_v0, tau_p0, - q_cfl.data()[eN_k]); + q_cfl(eN, k)); calculateSubgridError_tau(Ct_sge,Cd_sge, G,G_dd_G,tr_G, @@ -2773,7 +2783,7 @@ namespace proteus dmom_u_ham_grad_p[0], tau_v1, tau_p1, - q_cfl.data()[eN_k]); + q_cfl(eN, k)); tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0; tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0; @@ -2789,14 +2799,14 @@ namespace proteus subgridError_v, subgridError_w); // velocity used in adjoint (VMS or RBLES, with or without lagging the grid scale velocity) - dmom_adv_star[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*xt + useRBLES*subgridError_u); - dmom_adv_star[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v); + dmom_adv_star[0] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 0) - MOVING_DOMAIN*xt + useRBLES*subgridError_u); + dmom_adv_star[1] = inertial_term*dmom_u_acc_u*(q_velocity_sge(eN, k, 1) - MOVING_DOMAIN*yt + useRBLES*subgridError_v); - mom_u_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+0]); - mom_u_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+0]); + mom_u_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge(eN, k, 0)); + mom_u_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge(eN, k, 0)); - mom_v_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+1]); - mom_v_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+1]); + mom_v_adv[0] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge(eN, k, 1)); + mom_v_adv[1] += inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge(eN, k, 1)); // adjoint times the test functions for (int i=0;i 0) { //cek todo, this needs to be fixed for not exact @@ -2848,10 +2858,10 @@ namespace proteus { for (int I=0;I 0) if (element_active) { - q_mom_u_acc.data()[eN_k] = mom_u_acc; - q_mom_v_acc.data()[eN_k] = mom_v_acc; + q_mom_u_acc(eN, k) = mom_u_acc; + q_mom_v_acc(eN, k) = mom_v_acc; //subgrid error uses grid scale velocity - q_mass_adv.data()[eN_k_nSpace+0] = u; - q_mass_adv.data()[eN_k_nSpace+1] = v; + q_mass_adv(eN, k, 0) = u; + q_mass_adv(eN, k, 1) = v; } else//use the solid velocity { if (use_ball_as_particle) { - q_mom_u_acc.data()[eN_k] = particle_velocities.data()[eN_k_3d+0]; - q_mom_v_acc.data()[eN_k] = particle_velocities.data()[eN_k_3d+1]; - q_mass_adv.data()[eN_k_nSpace+0] = particle_velocities.data()[eN_k_3d+0]; - q_mass_adv.data()[eN_k_nSpace+1] = particle_velocities.data()[eN_k_3d+1]; + q_mom_u_acc(eN, k) = particle_velocities(eN, k, 0); + q_mom_v_acc(eN, k) = particle_velocities(eN, k, 1); + q_mass_adv(eN, k, 0) = particle_velocities(eN, k, 0); + q_mass_adv(eN, k, 1) = particle_velocities(eN, k, 1); } else { - q_mom_u_acc.data()[eN_k] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+0]; - q_mom_v_acc.data()[eN_k] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+1]; - q_mass_adv.data()[eN_k_nSpace+0] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+0]; - q_mass_adv.data()[eN_k_nSpace+1] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+1]; + q_mom_u_acc(eN, k) = particle_velocities(particle_index, eN, k, 0); + q_mom_v_acc(eN, k) = particle_velocities(particle_index, eN, k, 1); + q_mass_adv(eN, k, 0) = particle_velocities(particle_index, eN, k, 0); + q_mass_adv(eN, k, 1) = particle_velocities(particle_index, eN, k, 1); } } // @@ -3033,10 +3043,10 @@ namespace proteus v_L2 += v_e*v_e*H_s*dV*H_f; velocity_L2 += velocity_e*velocity_e*H_s*dV*H_f; p_dv += p*H_s*H_f*dV; - pa_dv += q_u_0.data()[eN_k]*H_s*H_f*dV; + pa_dv += q_u_0(eN, k)*H_s*H_f*dV; total_volume+=H_s*H_f*dV; total_surface_area+=D_s*H_f*dV; - if (phi_solid.data()[eN_k] >= 0.0) + if (phi_solid(eN, k) >= 0.0) { p_LI = fmax(p_LI, fabs(p_e)); u_LI = fmax(u_LI, fabs(u_e)); @@ -3048,18 +3058,18 @@ namespace proteus { register int i_nSpace=i*nSpace; elementResidual_mesh[i] += H_s*H_f*(ck.Reaction_weak(1.0,p_test_dV[i]) - - ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) - + ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last(eN, k)/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])); elementResidual_p[i] += H_s*H_f*(ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) + ck.Hamiltonian_weak(mass_ham, p_test_dV[i]) + DM*MOVING_DOMAIN*(ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) - - ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) - + ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last(eN, k)/dV) - ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])) + ck.Reaction_weak(mass_source,p_test_dV[i])); if (nDOF_test_element == nDOF_v_test_element) { elementResidual_p[i] += - H_s*H_f*(PRESSURE_PROJECTION_STABILIZATION * ck.pressureProjection_weak(mom_uu_diff_ten[1], p, p_element_avg, p_test_ref.data()[k*nDOF_test_element+i], dV) + + H_s*H_f*(PRESSURE_PROJECTION_STABILIZATION * ck.pressureProjection_weak(mom_uu_diff_ten[1], p, p_element_avg, p_test_ref(k, i), dV) + (1 - PRESSURE_PROJECTION_STABILIZATION) * ck.SubgridError(subgridError_u,Lstar_u_p[i]) + (1 - PRESSURE_PROJECTION_STABILIZATION) * ck.SubgridError(subgridError_v,Lstar_v_p[i])); } @@ -3084,7 +3094,7 @@ namespace proteus ck.Reaction_weak(mom_u_source+NONCONSERVATIVE_FORM*dmom_u_acc_u*u*div_mesh_velocity,vel_test_dV[i]) + ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) + MOMENTUM_SGE*VELOCITY_SGE*ck.SubgridError(subgridError_u,Lstar_u_u[i]) + - ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace])); + ck.NumericalDiffusion(q_numDiff_u_last(eN, k),grad_u,&vel_grad_test_dV[i_nSpace])); elementResidual_v[i] += H_s*H_f*(ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) + ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) + ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) + @@ -3092,7 +3102,7 @@ namespace proteus ck.Reaction_weak(mom_v_source+NONCONSERVATIVE_FORM*dmom_v_acc_v*v*div_mesh_velocity,vel_test_dV[i]) + ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) + MOMENTUM_SGE*VELOCITY_SGE*ck.SubgridError(subgridError_v,Lstar_v_v[i]) + - ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace])); + ck.NumericalDiffusion(q_numDiff_v_last(eN, k),grad_v,&vel_grad_test_dV[i_nSpace])); elementResidual_u[i] += H_s*H_f*MOMENTUM_SGE*PRESSURE_SGE*ck.SubgridError(subgridError_p,Lstar_p_u[i]); elementResidual_v[i] += H_s*H_f*MOMENTUM_SGE*PRESSURE_SGE*ck.SubgridError(subgridError_p,Lstar_p_v[i]); if (nParticles > 0)//solid boundary terms @@ -3106,12 +3116,12 @@ namespace proteus } }//i //estimate the numerical viscosity combining shock capturing and VMS/SUPG - numerical_viscosity.data()[eN_k] = q_numDiff_u_last.data()[eN_k] + MOMENTUM_SGE*VELOCITY_SGE*tau_v*(dmom_adv_star[0]*dmom_adv_star[0]+ + numerical_viscosity(eN, k) = q_numDiff_u_last(eN, k) + MOMENTUM_SGE*VELOCITY_SGE*tau_v*(dmom_adv_star[0]*dmom_adv_star[0]+ dmom_adv_star[1]*dmom_adv_star[1]); if (!isActiveElement[eN]) { - assert(std::fabs(gf_s.H(particle_eps,phi_solid.data()[eN_k])) == 0.0); - assert(std::fabs(gf_s.D(particle_eps,phi_solid.data()[eN_k])) == 0.0); + assert(std::fabs(gf_s.H(particle_eps,phi_solid(eN, k))) == 0.0); + assert(std::fabs(gf_s.D(particle_eps,phi_solid(eN, k))) == 0.0); } }//k }//fluid_phase @@ -3120,9 +3130,9 @@ namespace proteus { //compute indices and declare local storage register int eN_k = eN*nQuadraturePoints_element+k; - q_numDiff_u.data()[eN_k] = numDiffMax; - q_numDiff_v.data()[eN_k] = numDiffMax; - q_numDiff_w.data()[eN_k] = numDiffMax; + q_numDiff_u(eN, k) = numDiffMax; + q_numDiff_v(eN, k) = numDiffMax; + q_numDiff_w(eN, k) = numDiffMax; } #endif // @@ -3130,18 +3140,17 @@ namespace proteus // for(int i=0;i DWp_Dn_jump, DW_Dn_jump; - register double gamma_cutfem=ghost_penalty_constant,gamma_cutfem_p=ghost_penalty_constant,h_cutfem=elementBoundaryDiameter.data()[*it]; - int eN_nDOF_v_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_v_trial_element; + register double gamma_cutfem=ghost_penalty_constant,gamma_cutfem_p=ghost_penalty_constant,h_cutfem=elementBoundaryDiameter(*it); + int eN = elementBoundaryElementsArray.data()[(*it)*2+0]; + int eN_nDOF_v_trial_element = eN*nDOF_v_trial_element; //See Massing Schott Wall 2018 //cek todo modify for two-fluids: rho_0 != rho_1 double norm_v=0.0; for (int i_offset=1;i_offset(0.0); + auto elementResidual_p = detail::make_array(0.0); + auto elementResidual_u = detail::make_array(0.0); + auto elementResidual_v = detail::make_array(0.0); + register double eps_rho,eps_mu; const double* elementResidual_w(NULL); - for (int i=0;i 0) { - get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x_ext,y_ext,z_ext,ebqe_phi_s.data()[ebNE_kb]); + get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x_ext,y_ext,z_ext,ebqe_phi_s(ebNE, kb)); } //else ebqe_phi_s.data()[ebNE_kb] is computed in Prestep const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]); @@ -3653,19 +3654,19 @@ namespace proteus sigma, rho, nu, - elementDiameter.data()[eN], + elementDiameter(eN), smagorinskyConstant, turbulenceClosureModel, g.data(), useVF, - ebqe_vf_ext.data()[ebNE_kb], - ebqe_phi_ext.data()[ebNE_kb], - &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace], - ebqe_kappa_phi_ext.data()[ebNE_kb], + ebqe_vf_ext(ebNE, kb), + ebqe_phi_ext(ebNE, kb), + &ebqe_normal_phi_ext(ebNE, kb, 0), + ebqe_kappa_phi_ext(ebNE, kb), //VRANS porosity_ext, // - ebqe_phi_s.data()[ebNE_kb], + ebqe_phi_s(ebNE, kb), p_old, u_old, v_old, @@ -3683,8 +3684,8 @@ namespace proteus v_ext, w_ext, LAG_LES, - ebqe_eddy_viscosity.data()[ebNE_kb], - ebqe_eddy_viscosity_last.data()[ebNE_kb], + ebqe_eddy_viscosity(ebNE, kb), + ebqe_eddy_viscosity_last(ebNE, kb), mom_u_acc_ext, dmom_u_acc_u_ext, mom_v_acc_ext, @@ -3750,19 +3751,19 @@ namespace proteus sigma, rho, nu, - elementDiameter.data()[eN], + elementDiameter(eN), smagorinskyConstant, turbulenceClosureModel, g.data(), useVF, - bc_ebqe_vf_ext.data()[ebNE_kb], - bc_ebqe_phi_ext.data()[ebNE_kb], - &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace], - ebqe_kappa_phi_ext.data()[ebNE_kb], + bc_ebqe_vf_ext(ebNE, kb), + bc_ebqe_phi_ext(ebNE, kb), + &ebqe_normal_phi_ext(ebNE, kb, 0), + ebqe_kappa_phi_ext(ebNE, kb), //VRANS porosity_ext, // - ebqe_phi_s.data()[ebNE_kb], + ebqe_phi_s(ebNE, kb), p_old, u_old, v_old, @@ -3781,7 +3782,7 @@ namespace proteus bc_w_ext, LAG_LES, bc_eddy_viscosity_ext, - ebqe_eddy_viscosity_last.data()[ebNE_kb], + ebqe_eddy_viscosity_last(ebNE, kb), bc_mom_u_acc_ext, bc_dmom_u_acc_u_ext, bc_mom_v_acc_ext, @@ -3852,14 +3853,14 @@ namespace proteus rho_1, nu_1, useVF, - ebqe_vf_ext.data()[ebNE_kb], - ebqe_phi_ext.data()[ebNE_kb], + ebqe_vf_ext(ebNE, kb), + ebqe_phi_ext(ebNE, kb), porosity_ext, c_mu, //mwf hack - ebqe_turb_var_0.data()[ebNE_kb], - ebqe_turb_var_1.data()[ebNE_kb], + ebqe_turb_var_0(ebNE, kb), + ebqe_turb_var_1(ebNE, kb), turb_var_grad_0_dummy, //not needed - ebqe_eddy_viscosity.data()[ebNE_kb], + ebqe_eddy_viscosity(ebNE, kb), mom_uu_diff_ten_ext, mom_vv_diff_ten_ext, mom_ww_diff_ten_ext, @@ -3882,12 +3883,12 @@ namespace proteus rho_1, nu_1, useVF, - bc_ebqe_vf_ext.data()[ebNE_kb], - bc_ebqe_phi_ext.data()[ebNE_kb], + bc_ebqe_vf_ext(ebNE, kb), + bc_ebqe_phi_ext(ebNE, kb), porosity_ext, c_mu, //mwf hack - ebqe_turb_var_0.data()[ebNE_kb], - ebqe_turb_var_1.data()[ebNE_kb], + ebqe_turb_var_0(ebNE, kb), + ebqe_turb_var_1(ebNE, kb), turb_var_grad_0_dummy, //not needed bc_eddy_viscosity_ext, bc_mom_uu_diff_ten_ext, @@ -3950,16 +3951,16 @@ namespace proteus //calculate the numerical fluxes // ck.calculateGScale(G,normal,h_penalty); - penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb]; + penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext(ebNE, kb); exteriorNumericalAdvectiveFlux(NONCONSERVATIVE_FORM, - isDOFBoundary_p.data()[ebNE_kb], - isDOFBoundary_u.data()[ebNE_kb], - isDOFBoundary_v.data()[ebNE_kb], - isDOFBoundary_w.data()[ebNE_kb], - isAdvectiveFluxBoundary_p.data()[ebNE_kb], - isAdvectiveFluxBoundary_u.data()[ebNE_kb], - isAdvectiveFluxBoundary_v.data()[ebNE_kb], - isAdvectiveFluxBoundary_w.data()[ebNE_kb], + isDOFBoundary_p(ebNE, kb), + isDOFBoundary_u(ebNE, kb), + isDOFBoundary_v(ebNE, kb), + isDOFBoundary_w(ebNE, kb), + isAdvectiveFluxBoundary_p(ebNE, kb), + isAdvectiveFluxBoundary_u(ebNE, kb), + isAdvectiveFluxBoundary_v(ebNE, kb), + isAdvectiveFluxBoundary_w(ebNE, kb), dmom_u_ham_grad_p_ext[0],//=1/rho, bc_dmom_u_ham_grad_p_ext[0],//=1/bc_rho, normal, @@ -3970,10 +3971,10 @@ namespace proteus bc_mom_u_adv_ext, bc_mom_v_adv_ext, bc_mom_w_adv_ext, - ebqe_bc_flux_mass_ext.data()[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),//BC is relative mass flux - ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb], - ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb], - ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb], + ebqe_bc_flux_mass_ext(ebNE, kb)+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),//BC is relative mass flux + ebqe_bc_flux_mom_u_adv_ext(ebNE, kb), + ebqe_bc_flux_mom_v_adv_ext(ebNE, kb), + ebqe_bc_flux_mom_w_adv_ext(ebNE, kb), p_ext, u_ext, v_ext, @@ -4001,30 +4002,30 @@ namespace proteus flux_mom_u_adv_ext, flux_mom_v_adv_ext, flux_mom_w_adv_ext, - &ebqe_velocity.data()[ebNE_kb_nSpace]); + &ebqe_velocity(ebNE, kb, 0)); for (int I=0;I nDOF_pressure) - nDOF_pressure=p_l2g.data()[eN_i]; + auto tmp = p_l2g(eN, i); + if (tmp > nDOF_pressure) + nDOF_pressure=tmp; } } nDOF_pressure +=1; assert(p_dof.shape(0) == nDOF_pressure); //std::cout<<"nDOF_pressure "< p_LI)