18 #define TURB_FORCE_FAC 0.0
19 #define CUT_CELL_INTEGRATION 0
20 double sgn(
double val) {
21 return double((0.0 < val) - (val < 0.0));
37 #define CELL_BASED_EV_COEFF 1
38 #define POWER_SMOOTHNESS_INDICATOR 2
39 #define EPS_FOR_GAMMA_INDICATOR 1E-10
40 #define C_FOR_GAMMA_INDICATOR 0.25 // increase gamma to make the indicator more agressive (less dissipative)
41 #define USE_GAMMA_INDICATOR 0
42 #define ANISOTROPIC_DIFFUSION 0
50 double detT = (r1[1] - r2[1])*(r0[0] - r2[0]) + (r2[0] - r1[0])*(r0[1] - r2[1]);
51 lambda[0] = ((r1[1] - r2[1])*(
r[0] - r2[0]) + (r2[0] - r1[0])*(
r[1] - r2[1]))/detT;
52 lambda[1] = ((r2[1] - r0[1])*(
r[0] - r2[0]) + (r0[0] - r2[0])*(
r[1] - r2[1]))/detT;
53 lambda[2] = 1.0 - lambda[0] - lambda[1];
58 template<
int nSpace,
int nP,
int nQ,
int nEBQ>
98 template<
class CompKernelType,
100 int nQuadraturePoints_element,
101 int nDOF_mesh_trial_element,
102 int nDOF_trial_element,
103 int nDOF_test_element,
104 int nQuadraturePoints_elementBoundary>
133 M_PI/6., 0.05, 1.00),
168 double mu_fr_limiter)
189 inline double Dot(
const double vec1[nSpace],
190 const double vec2[nSpace])
193 for (
int I=0; I<nSpace; I++)
194 dot += vec1[I]*vec2[I];
199 const double vel_grad[nSpace],
200 double vel_tgrad[nSpace])
202 double normal_dot_vel_grad =
Dot(normal,vel_grad);
203 for (
int I=0; I<nSpace; I++)
204 vel_tgrad[I] = vel_grad[I] - normal_dot_vel_grad*normal[I];
217 const double smagorinskyConstant,
218 const int turbulenceClosureModel,
219 const double g[nSpace],
223 const double n[nSpace],
224 const double distance_to_omega_solid,
226 const double porosity,
228 const double grad_p[nSpace],
229 const double grad_u[nSpace],
230 const double grad_v[nSpace],
231 const double grad_w[nSpace],
238 double& eddy_viscosity,
240 double& dmom_u_acc_u,
242 double& dmom_v_acc_v,
244 double& dmom_w_acc_w,
245 double mass_adv[nSpace],
246 double dmass_adv_u[nSpace],
247 double dmass_adv_v[nSpace],
248 double dmass_adv_w[nSpace],
249 double mom_u_adv[nSpace],
250 double dmom_u_adv_u[nSpace],
251 double dmom_u_adv_v[nSpace],
252 double dmom_u_adv_w[nSpace],
253 double mom_v_adv[nSpace],
254 double dmom_v_adv_u[nSpace],
255 double dmom_v_adv_v[nSpace],
256 double dmom_v_adv_w[nSpace],
257 double mom_w_adv[nSpace],
258 double dmom_w_adv_u[nSpace],
259 double dmom_w_adv_v[nSpace],
260 double dmom_w_adv_w[nSpace],
261 double mom_uu_diff_ten[nSpace],
262 double mom_vv_diff_ten[nSpace],
263 double mom_ww_diff_ten[nSpace],
264 double mom_uv_diff_ten[1],
265 double mom_uw_diff_ten[1],
266 double mom_vu_diff_ten[1],
267 double mom_vw_diff_ten[1],
268 double mom_wu_diff_ten[1],
269 double mom_wv_diff_ten[1],
270 double& mom_u_source,
271 double& mom_v_source,
272 double& mom_w_source,
274 double dmom_u_ham_grad_p[nSpace],
275 double dmom_u_ham_grad_u[nSpace],
277 double dmom_v_ham_grad_p[nSpace],
278 double dmom_v_ham_grad_v[nSpace],
280 double dmom_w_ham_grad_p[nSpace],
281 double dmom_w_ham_grad_w[nSpace],
284 int KILL_PRESSURE_TERM,
285 int MULTIPLY_EXTERNAL_FORCE_BY_DENSITY,
289 int MATERIAL_PARAMETERS_AS_FUNCTION,
290 double density_as_function,
291 double dynamic_viscosity_as_function,
293 double x,
double y,
double z,
294 int use_ball_as_particle,
297 double* ball_velocity,
298 double* ball_angular_velocity,
300 int INT_BY_PARTS_PRESSURE)
302 double rho,nu,mu,H_rho,ImH_rho,d_rho,H_mu,ImH_mu,d_mu,norm_n,nu_t0=0.0,nu_t1=0.0,nu_t;
303 H_rho = (1.0-useVF)*
gf.H(eps_rho,
phi) + useVF*fmin(1.0,fmax(0.0,vf));
304 ImH_rho = (1.0-useVF)*
gf.ImH(eps_rho,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
305 d_rho = (1.0-useVF)*
gf.D(eps_rho,
phi);
306 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi) + useVF*fmin(1.0,fmax(0.0,vf));
307 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi) + useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
308 d_mu = (1.0-useVF)*
gf.D(eps_mu,
phi);
311 switch (turbulenceClosureModel)
316 norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] +
317 0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0])));
319 nu_t0 = smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
320 nu_t1 = smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
324 double re_0,cs_0=0.0,re_1,cs_1=0.0;
325 norm_S = sqrt(2.0*(grad_u[0]*grad_u[0] + grad_v[1]*grad_v[1] +
326 0.5*(grad_u[1]+grad_v[0])*(grad_u[1]+grad_v[0])));
327 re_0 = h_e*h_e*norm_S/
nu_0;
329 cs_0=0.027*pow(10.0,-3.23*pow(re_0,-0.92));
330 nu_t0 = cs_0*h_e*h_e*norm_S;
331 re_1 = h_e*h_e*norm_S/
nu_1;
333 cs_1=0.027*pow(10.0,-3.23*pow(re_1,-0.92));
334 nu_t1 = cs_1*h_e*h_e*norm_S;
338 if (MATERIAL_PARAMETERS_AS_FUNCTION==0)
341 nu_t= nu_t0*ImH_mu+nu_t1*H_mu;
348 rho = density_as_function;
350 mu = dynamic_viscosity_as_function;
357 eddy_viscosity = nu_t*rho;
361 double phi_s_effect = (distance_to_omega_solid > 0.0) ? 1.0 : 1e-10;
367 dmom_u_acc_u=phi_s_effect * rho*porosity;
371 dmom_v_acc_v=phi_s_effect * rho*porosity;
378 mass_adv[0]=phi_s_effect * porosity*
u;
379 mass_adv[1]=phi_s_effect * porosity*
v;
382 dmass_adv_u[0]=phi_s_effect * porosity;
387 dmass_adv_v[1]=phi_s_effect * porosity;
447 mom_uu_diff_ten[0] = phi_s_effect * porosity*2.0*mu;
448 mom_uu_diff_ten[1] = phi_s_effect * porosity*mu;
451 mom_uv_diff_ten[0]=phi_s_effect * porosity*mu;
456 mom_vv_diff_ten[0] = phi_s_effect * porosity*mu;
457 mom_vv_diff_ten[1] = phi_s_effect * porosity*2.0*mu;
460 mom_vu_diff_ten[0]=phi_s_effect * porosity*mu;
474 norm_n = sqrt(
n[0]*
n[0]+
n[1]*
n[1]);
475 mom_u_source = -phi_s_effect * porosity*rho*g[0];
476 mom_v_source = -phi_s_effect * porosity*rho*g[1];
480 mom_u_source -= (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcex;
481 mom_v_source -= (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcey;
485 double aux_pressure = (KILL_PRESSURE_TERM==1 ? 0. : 1.)*(INT_BY_PARTS_PRESSURE==1 ? 0. : 1.);
486 mom_u_ham = phi_s_effect * porosity*grad_p[0]*aux_pressure;
487 dmom_u_ham_grad_p[0]=phi_s_effect * porosity*aux_pressure;
488 dmom_u_ham_grad_p[1]=0.0;
492 mom_v_ham = phi_s_effect * porosity*grad_p[1]*aux_pressure;
493 dmom_v_ham_grad_p[0]=0.0;
494 dmom_v_ham_grad_p[1]=phi_s_effect * porosity*aux_pressure;
504 mom_u_ham += phi_s_effect * porosity*rho*(uStar*grad_u[0]+vStar*grad_u[1]);
505 dmom_u_ham_grad_u[0]=phi_s_effect * porosity*rho*uStar;
506 dmom_u_ham_grad_u[1]=phi_s_effect * porosity*rho*vStar;
510 mom_v_ham += phi_s_effect * porosity*rho*(uStar*grad_v[0]+vStar*grad_v[1]);
511 dmom_v_ham_grad_v[0]=phi_s_effect * porosity*rho*uStar;
512 dmom_v_ham_grad_v[1]=phi_s_effect * porosity*rho*vStar;
530 const double eps_rho,
551 const double uStar_s,
552 const double vStar_s,
553 const double wStar_s,
554 double& mom_u_source,
555 double& mom_v_source,
556 double& mom_w_source,
557 double dmom_u_source[nSpace],
558 double dmom_v_source[nSpace],
559 double dmom_w_source[nSpace],
564 double rho, mu,nu,H_mu,ImH_mu,uc,duc_du,duc_dv,duc_dw,viscosity,H_s;
565 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi)+useVF*fmin(1.0,fmax(0.0,vf));
566 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi)+useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
572 duc_du =
u/(uc+1.0e-12);
573 duc_dv =
v/(uc+1.0e-12);
574 duc_dw =
w/(uc+1.0e-12);
575 double fluid_velocity[2]={uStar,vStar}, solid_velocity[2]={uStar_s,vStar_s};
582 double beta2 = 156976.4;
585 (1.0 - phi_s)*(1.0-
DRAG_FAC)*beta2*(
u-u_s);
587 (1.0 - phi_s)*(1.0-
DRAG_FAC)*beta2*(
v-v_s);
591 dmom_u_source[0] = (1.0 - phi_s) * new_beta + (1.0 - phi_s)*(1.0-
DRAG_FAC)*beta2;
592 dmom_u_source[1] = 0.0;
595 dmom_v_source[0] = 0.0;
596 dmom_v_source[1] = (1.0 - phi_s) * new_beta + (1.0 - phi_s)*(1.0-
DRAG_FAC)*beta2;
599 dmom_w_source[0] = 0.0;
600 dmom_w_source[1] = 0.0;
605 const double particle_nitsche,
607 const int nParticles,
609 double *particle_signed_distances,
610 double *particle_signed_distance_normals,
611 double *particle_velocities,
612 double *particle_centroids,
613 int use_ball_as_particle,
616 double* ball_velocity,
617 double* ball_angular_velocity,
618 const double porosity,
619 const double penalty,
622 const double eps_rho,
642 const double grad_u[nSpace],
643 const double grad_v[nSpace],
644 const double grad_w[nSpace],
645 double &mom_u_source,
646 double &mom_v_source,
647 double &mom_w_source,
648 double dmom_u_source[nSpace],
649 double dmom_v_source[nSpace],
650 double dmom_w_source[nSpace],
651 double mom_u_adv[nSpace],
652 double mom_v_adv[nSpace],
653 double mom_w_adv[nSpace],
654 double dmom_u_adv_u[nSpace],
655 double dmom_v_adv_v[nSpace],
656 double dmom_w_adv_w[nSpace],
658 double dmom_u_ham_grad_u[nSpace],
660 double dmom_v_ham_grad_v[nSpace],
662 double dmom_w_ham_grad_w[nSpace],
663 double *particle_netForces,
664 double *particle_netMoments,
665 double *particle_surfaceArea)
667 double C, rho, mu, nu, H_mu, ImH_mu, uc, duc_du, duc_dv, duc_dw, H_s, D_s, phi_s, u_s, v_s, w_s;
668 double force_x, force_y, r_x, r_y, force_p_x, force_p_y, force_stress_x, force_stress_y;
669 double phi_s_normal[2]={0.0};
670 double fluid_outward_normal[2];
673 H_mu = (1.0 - useVF) *
gf.H(eps_mu,
phi) + useVF * fmin(1.0, fmax(0.0, vf));
674 ImH_mu = (1.0 - useVF) *
gf.ImH(eps_mu,
phi) + useVF * (1.0-fmin(1.0, fmax(0.0, vf)));
679 for (
int i = 0; i < nParticles; i++)
681 double* vel_pointer= &particle_velocities[i*sd_offset*nSpace];
682 if(use_ball_as_particle==1)
687 ball_velocity,ball_angular_velocity,
690 center[0] = ball_center[3*i+0];
691 center[1] = ball_center[3*i+1];
697 phi_s = particle_signed_distances[i * sd_offset];
698 phi_s_normal[0] = particle_signed_distance_normals[i * sd_offset * 3 + 0];
699 phi_s_normal[1] = particle_signed_distance_normals[i * sd_offset * 3 + 1];
700 vel[0] = particle_velocities[i * sd_offset * 3 + 0];
701 vel[1] = particle_velocities[i * sd_offset * 3 + 1];
702 u_s = vel_pointer[0];
703 v_s = vel_pointer[1];
704 center[0] = particle_centroids[3*i+0];
705 center[1] = particle_centroids[3*i+1];
708 fluid_outward_normal[0] = -phi_s_normal[0];
709 fluid_outward_normal[1] = -phi_s_normal[1];
712 H_s =
gf_s.H(eps_s, phi_s);
713 D_s =
gf_s.D(eps_s, phi_s);
714 double rel_vel_norm = sqrt((uStar - u_s) * (uStar - u_s) +
715 (vStar - v_s) * (vStar - v_s) +
716 (wStar - w_s) * (wStar - w_s));
720 double C_surf = mu * penalty;
721 double C_vol = (alpha + beta * rel_vel_norm);
723 C = (D_s * C_surf +
gf_s.ImH(eps_s, phi_s) * C_vol);
724 force_x = dV * D_s * (p * fluid_outward_normal[0]
725 - porosity*mu*(fluid_outward_normal[0] * 2* grad_u[0] +
726 fluid_outward_normal[1] * (grad_u[1]+grad_v[0]))
729 force_y = dV * D_s * (p * fluid_outward_normal[1]
730 - porosity*mu * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) +
731 fluid_outward_normal[1] * 2* grad_v[1])
734 force_p_x = dV * D_s * p * fluid_outward_normal[0];
735 force_p_y = dV * D_s * p * fluid_outward_normal[1];
736 force_stress_x = dV * D_s * (-porosity*mu * (fluid_outward_normal[0] * 2* grad_u[0] +
737 fluid_outward_normal[1] * (grad_u[1]+grad_v[0]))
740 force_stress_y = dV * D_s * (-porosity*mu * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) +
741 fluid_outward_normal[1] * 2* grad_v[1])
750 particle_surfaceArea[i] += dV * D_s;
751 particle_netForces[i * 3 + 0] += force_x;
752 particle_netForces[i * 3 + 1] += force_y;
753 particle_netForces[(i+ nParticles)*3+0]+= force_p_x;
754 particle_netForces[(i+2*nParticles)*3+0]+= force_stress_x;
755 particle_netForces[(i+ nParticles)*3+1]+= force_p_y;
756 particle_netForces[(i+2*nParticles)*3+1]+= force_stress_y;
757 particle_netMoments[i * 3 + 2] += (r_x * force_y - r_y * force_x);
761 mom_u_source += C * (
u - u_s);
762 mom_v_source += C * (
v - v_s);
764 dmom_u_source[0] += C;
765 dmom_v_source[1] += C;
768 mom_u_ham -= D_s * porosity * mu * (fluid_outward_normal[0] * grad_u[0] + fluid_outward_normal[1] * grad_u[1]);
769 dmom_u_ham_grad_u[0] -= D_s * porosity * mu * fluid_outward_normal[0];
770 dmom_u_ham_grad_u[1] -= D_s * porosity * mu * fluid_outward_normal[1];
772 mom_v_ham -= D_s * porosity * mu * (fluid_outward_normal[0] * grad_v[0] + fluid_outward_normal[1] * grad_v[1]);
773 dmom_v_ham_grad_v[0] -= D_s * porosity * mu * fluid_outward_normal[0];
774 dmom_v_ham_grad_v[1] -= D_s * porosity * mu * fluid_outward_normal[1];
776 mom_u_adv[0] += D_s * porosity * mu * fluid_outward_normal[0] * (
u - u_s);
777 mom_u_adv[1] += D_s * porosity * mu * fluid_outward_normal[1] * (
u - u_s);
778 dmom_u_adv_u[0] += D_s * porosity * mu * fluid_outward_normal[0];
779 dmom_u_adv_u[1] += D_s * porosity * mu * fluid_outward_normal[1];
781 mom_v_adv[0] += D_s * porosity * mu * fluid_outward_normal[0] * (
v - v_s);
782 mom_v_adv[1] += D_s * porosity * mu * fluid_outward_normal[1] * (
v - v_s);
783 dmom_v_adv_v[0] += D_s * porosity * mu * fluid_outward_normal[0];
784 dmom_v_adv_v[1] += D_s * porosity * mu * fluid_outward_normal[1];
789 const int nParticles,
791 double *particle_signed_distances,
792 double *particle_signed_distance_normals,
793 double *particle_velocities,
794 double *particle_centroids,
795 int use_ball_as_particle,
798 double* ball_velocity,
799 double* ball_angular_velocity,
800 const double penalty,
803 const double eps_rho,
823 const double grad_u[nSpace],
824 const double grad_v[nSpace],
825 const double grad_w[nSpace],
826 double* particle_netForces,
827 double* particle_netMoments)
829 double C, rho, mu, nu, H_mu, ImH_mu, uc, duc_du, duc_dv, duc_dw, H_s, D_s, phi_s, u_s, v_s, w_s, force_x, force_y, r_x, r_y;
830 double phi_s_normal[2];
831 double fluid_outward_normal[2];
834 H_mu = (1.0 - useVF) *
gf.H(eps_mu,
phi) + useVF * fmin(1.0, fmax(0.0, vf));
835 ImH_mu = (1.0 - useVF) *
gf.ImH(eps_mu,
phi) + useVF * (1.0-fmin(1.0, fmax(0.0, vf)));
840 for (
int i = 0; i < nParticles; i++)
842 if(use_ball_as_particle==1)
847 ball_velocity,ball_angular_velocity,
850 center[0] = ball_center[3*i+0];
851 center[1] = ball_center[3*i+1];
855 phi_s = particle_signed_distances[i * sd_offset];
856 phi_s_normal[0] = particle_signed_distance_normals[i * sd_offset * 3 + 0];
857 phi_s_normal[1] = particle_signed_distance_normals[i * sd_offset * 3 + 1];
858 vel[0] = particle_velocities[i * sd_offset * 3 + 0];
859 vel[1] = particle_velocities[i * sd_offset * 3 + 1];
860 center[0] = particle_centroids[3*i+0];
861 center[1] = particle_centroids[3*i+1];
864 fluid_outward_normal[0] = -phi_s_normal[0];
865 fluid_outward_normal[1] = -phi_s_normal[1];
869 H_s =
gf_s.H(eps_s, phi_s);
870 D_s =
gf_s.D(eps_s, phi_s);
871 double rel_vel_norm = sqrt((uStar - u_s) * (uStar - u_s) +(vStar - v_s) * (vStar - v_s) + (wStar - w_s) * (wStar - w_s));
872 double C_surf = (phi_s > 0.0) ? 0.0 : nu * penalty;
873 double C_vol = (phi_s > 0.0) ? 0.0 : (alpha + beta * rel_vel_norm);
875 C = (D_s * C_surf +
gf_s.ImH(eps_s, phi_s) * C_vol);
876 force_x = dV * D_s * (p * fluid_outward_normal[0]
877 -mu * (fluid_outward_normal[0] * 2* grad_u[0] + fluid_outward_normal[1] * (grad_u[1]+grad_v[0]))
881 force_y = dV * D_s * (p * fluid_outward_normal[1]
882 -mu * (fluid_outward_normal[0] * (grad_u[1]+grad_v[0]) + fluid_outward_normal[1] * 2* grad_v[1])
893 particle_netForces[i * 3 + 0] += force_x;
894 particle_netForces[i * 3 + 1] += force_y;
895 particle_netMoments[i * 3 + 2] += (r_x * force_y - r_y * force_x);
901 const double& elementDiameter,
903 const double df[nSpace],
906 double h,density,nrm_df=0.0;
907 h = hFactor*elementDiameter;
909 for(
int I=0;I<nSpace;I++)
911 nrm_df = sqrt(nrm_df);
912 if (density > 1.0e-8)
913 cfl = nrm_df/(h*density);
920 const double eps_rho,
929 const double porosity,
930 const double eddy_visc_coef_0,
931 const double turb_var_0,
932 const double turb_var_1,
933 const double turb_grad_0[nSpace],
934 double &eddy_viscosity,
935 double mom_uu_diff_ten[nSpace],
936 double mom_vv_diff_ten[nSpace],
937 double mom_ww_diff_ten[nSpace],
938 double mom_uv_diff_ten[1],
939 double mom_uw_diff_ten[1],
940 double mom_vu_diff_ten[1],
941 double mom_vw_diff_ten[1],
942 double mom_wu_diff_ten[1],
943 double mom_wv_diff_ten[1],
944 double &mom_u_source,
945 double &mom_v_source,
946 double &mom_w_source)
954 assert (turbulenceClosureModel >=3);
955 double rho,nu,H_mu,ImH_mu,nu_t=0.0,nu_t_keps =0.0, nu_t_komega=0.0;
956 double isKEpsilon = 1.0;
957 if (turbulenceClosureModel == 4)
959 H_mu = (1.0-useVF)*
gf.H(eps_mu,
phi)+useVF*fmin(1.0,fmax(0.0,vf));
960 ImH_mu = (1.0-useVF)*
gf.ImH(eps_mu,
phi)+useVF*(1.0-fmin(1.0,fmax(0.0,vf)));
964 const double twoThirds = 2.0/3.0;
const double div_zero = 1.0e-2*fmin(
nu_0,
nu_1);
965 mom_u_source += twoThirds*turb_grad_0[0];
966 mom_v_source += twoThirds*turb_grad_0[1];
971 nu_t_keps = eddy_visc_coef_0*turb_var_0*turb_var_0/(fabs(turb_var_1) + div_zero);
973 nu_t_komega = turb_var_0/(fabs(turb_var_1) + div_zero);
975 nu_t = isKEpsilon*nu_t_keps + (1.0-isKEpsilon)*nu_t_komega;
982 nu_t = fmax(nu_t,1.0e-4*nu);
984 nu_t = fmin(nu_t,1.0e6*nu);
986 eddy_viscosity = nu_t*rho;
988 mom_uu_diff_ten[0] += porosity*2.0*eddy_viscosity;
989 mom_uu_diff_ten[1] += porosity*eddy_viscosity;
992 mom_uv_diff_ten[0]+=porosity*eddy_viscosity;
997 mom_vv_diff_ten[0] += porosity*eddy_viscosity;
998 mom_vv_diff_ten[1] += porosity*2.0*eddy_viscosity;
1001 mom_vu_diff_ten[0]+=porosity*eddy_viscosity;
1016 const double &elementDiameter,
1019 const double df[nSpace],
1026 double h, oneByAbsdt, density, viscosity, nrm_df;
1027 h = hFactor * elementDiameter;
1031 for (
int I = 0; I < nSpace; I++)
1032 nrm_df +=
df[I] *
df[I];
1033 nrm_df = sqrt(nrm_df);
1034 if (density > 1.0e-8)
1035 cfl = nrm_df/(h*density);
1038 oneByAbsdt = fabs(dmt);
1039 tau_v = 1.0/(4.0*viscosity/(h*h) + 2.0*nrm_df/h + oneByAbsdt);
1040 tau_p = (4.0*viscosity + 2.0*nrm_df*h + oneByAbsdt*h*h)/pfac;
1044 const double &Cd_sge,
1045 const double G[nSpace * nSpace],
1046 const double &G_dd_G,
1049 const double Ai[nSpace],
1056 double v_d_Gv = 0.0;
1057 for (
int I = 0; I < nSpace; I++)
1058 for (
int J = 0; J < nSpace; J++)
1059 v_d_Gv += Ai[I] * G[I * nSpace + J] * Ai[J];
1060 tau_v = 1.0 / sqrt(Ct_sge * A0 * A0 + v_d_Gv + Cd_sge * Kij * Kij * G_dd_G + 1.0e-12);
1061 tau_p = 1.0 / (pfac * tr_G * tau_v);
1065 const double &tau_v,
1066 const double &pdeResidualP,
1067 const double &pdeResidualU,
1068 const double &pdeResidualV,
1069 const double &pdeResidualW,
1070 double &subgridErrorP,
1071 double &subgridErrorU,
1072 double &subgridErrorV,
1073 double &subgridErrorW)
1076 subgridErrorP = -tau_p * pdeResidualP;
1078 subgridErrorU = -tau_v * pdeResidualU;
1079 subgridErrorV = -tau_v * pdeResidualV;
1084 const double &tau_v,
1085 const double dpdeResidualP_du[nDOF_trial_element],
1086 const double dpdeResidualP_dv[nDOF_trial_element],
1087 const double dpdeResidualP_dw[nDOF_trial_element],
1088 const double dpdeResidualU_dp[nDOF_trial_element],
1089 const double dpdeResidualU_du[nDOF_trial_element],
1090 const double dpdeResidualV_dp[nDOF_trial_element],
1091 const double dpdeResidualV_dv[nDOF_trial_element],
1092 const double dpdeResidualW_dp[nDOF_trial_element],
1093 const double dpdeResidualW_dw[nDOF_trial_element],
1094 double dsubgridErrorP_du[nDOF_trial_element],
1095 double dsubgridErrorP_dv[nDOF_trial_element],
1096 double dsubgridErrorP_dw[nDOF_trial_element],
1097 double dsubgridErrorU_dp[nDOF_trial_element],
1098 double dsubgridErrorU_du[nDOF_trial_element],
1099 double dsubgridErrorV_dp[nDOF_trial_element],
1100 double dsubgridErrorV_dv[nDOF_trial_element],
1101 double dsubgridErrorW_dp[nDOF_trial_element],
1102 double dsubgridErrorW_dw[nDOF_trial_element])
1104 for (
int j = 0; j < nDOF_trial_element; j++)
1107 dsubgridErrorP_du[j] = -tau_p * dpdeResidualP_du[j];
1108 dsubgridErrorP_dv[j] = -tau_p * dpdeResidualP_dv[j];
1112 dsubgridErrorU_dp[j] = -tau_v * dpdeResidualU_dp[j];
1113 dsubgridErrorU_du[j] = -tau_v * dpdeResidualU_du[j];
1115 dsubgridErrorV_dp[j] = -tau_v * dpdeResidualV_dp[j];
1116 dsubgridErrorV_dv[j] = -tau_v * dpdeResidualV_dv[j];
1125 const int& isDOFBoundary_u,
1126 const int& isDOFBoundary_v,
1127 const int& isDOFBoundary_w,
1128 const int& isFluxBoundary_p,
1129 const int& isFluxBoundary_u,
1130 const int& isFluxBoundary_v,
1131 const int& isFluxBoundary_w,
1132 const double& oneByRho,
1133 const double& bc_oneByRho,
1134 const double n[nSpace],
1135 const double& porosity,
1140 const double bc_f_mass[nSpace],
1141 const double bc_f_umom[nSpace],
1142 const double bc_f_vmom[nSpace],
1143 const double bc_f_wmom[nSpace],
1144 const double& bc_flux_mass,
1145 const double& bc_flux_umom,
1146 const double& bc_flux_vmom,
1147 const double& bc_flux_wmom,
1152 const double f_mass[nSpace],
1153 const double f_umom[nSpace],
1154 const double f_vmom[nSpace],
1155 const double f_wmom[nSpace],
1156 const double df_mass_du[nSpace],
1157 const double df_mass_dv[nSpace],
1158 const double df_mass_dw[nSpace],
1159 const double df_umom_dp[nSpace],
1160 const double df_umom_du[nSpace],
1161 const double df_umom_dv[nSpace],
1162 const double df_umom_dw[nSpace],
1163 const double df_vmom_dp[nSpace],
1164 const double df_vmom_du[nSpace],
1165 const double df_vmom_dv[nSpace],
1166 const double df_vmom_dw[nSpace],
1167 const double df_wmom_dp[nSpace],
1168 const double df_wmom_du[nSpace],
1169 const double df_wmom_dv[nSpace],
1170 const double df_wmom_dw[nSpace],
1175 double* velocity_star,
1178 double flowSpeedNormal;
1183 flowSpeedNormal=porosity*(
n[0]*velocity_star[0] +
1184 n[1]*velocity_star[1]);
1188 if (isDOFBoundary_u != 1)
1190 flux_mass +=
n[0]*f_mass[0];
1191 if (flowSpeedNormal < 0.0)
1193 flux_umom+=flowSpeedNormal*(0.0 -
u);
1198 flux_mass +=
n[0]*f_mass[0];
1199 if (flowSpeedNormal < 0.0)
1201 flux_umom+=flowSpeedNormal*(bc_u -
u);
1205 if (isDOFBoundary_v != 1)
1207 flux_mass+=
n[1]*f_mass[1];
1208 if (flowSpeedNormal < 0.0)
1210 flux_vmom+=flowSpeedNormal*(0.0 -
v);
1215 flux_mass+=
n[1]*f_mass[1];
1216 if (flowSpeedNormal < 0.0)
1218 flux_vmom+=flowSpeedNormal*(bc_v -
v);
1244 if (isFluxBoundary_u == 1)
1246 flux_umom = bc_flux_umom;
1247 velocity[0] = bc_flux_umom/porosity;
1249 if (isFluxBoundary_v == 1)
1251 flux_vmom = bc_flux_vmom;
1252 velocity[1] = bc_flux_umom/porosity;
1262 const int& isDOFBoundary_u,
1263 const int& isDOFBoundary_v,
1264 const int& isDOFBoundary_w,
1265 const int& isFluxBoundary_p,
1266 const int& isFluxBoundary_u,
1267 const int& isFluxBoundary_v,
1268 const int& isFluxBoundary_w,
1269 const double& oneByRho,
1270 const double n[nSpace],
1271 const double& porosity,
1276 const double bc_f_mass[nSpace],
1277 const double bc_f_umom[nSpace],
1278 const double bc_f_vmom[nSpace],
1279 const double bc_f_wmom[nSpace],
1280 const double& bc_flux_mass,
1281 const double& bc_flux_umom,
1282 const double& bc_flux_vmom,
1283 const double& bc_flux_wmom,
1288 const double f_mass[nSpace],
1289 const double f_umom[nSpace],
1290 const double f_vmom[nSpace],
1291 const double f_wmom[nSpace],
1292 const double df_mass_du[nSpace],
1293 const double df_mass_dv[nSpace],
1294 const double df_mass_dw[nSpace],
1295 const double df_umom_dp[nSpace],
1296 const double df_umom_du[nSpace],
1297 const double df_umom_dv[nSpace],
1298 const double df_umom_dw[nSpace],
1299 const double df_vmom_dp[nSpace],
1300 const double df_vmom_du[nSpace],
1301 const double df_vmom_dv[nSpace],
1302 const double df_vmom_dw[nSpace],
1303 const double df_wmom_dp[nSpace],
1304 const double df_wmom_du[nSpace],
1305 const double df_wmom_dv[nSpace],
1306 const double df_wmom_dw[nSpace],
1307 double& dflux_mass_du,
1308 double& dflux_mass_dv,
1309 double& dflux_mass_dw,
1310 double& dflux_umom_dp,
1311 double& dflux_umom_du,
1312 double& dflux_umom_dv,
1313 double& dflux_umom_dw,
1314 double& dflux_vmom_dp,
1315 double& dflux_vmom_du,
1316 double& dflux_vmom_dv,
1317 double& dflux_vmom_dw,
1318 double& dflux_wmom_dp,
1319 double& dflux_wmom_du,
1320 double& dflux_wmom_dv,
1321 double& dflux_wmom_dw,
1322 double* velocity_star)
1324 double flowSpeedNormal;
1325 dflux_mass_du = 0.0;
1326 dflux_mass_dv = 0.0;
1329 dflux_umom_dp = 0.0;
1330 dflux_umom_du = 0.0;
1331 dflux_umom_dv = 0.0;
1334 dflux_vmom_dp = 0.0;
1335 dflux_vmom_du = 0.0;
1336 dflux_vmom_dv = 0.0;
1339 dflux_wmom_dp = 0.0;
1340 dflux_wmom_du = 0.0;
1341 dflux_wmom_dv = 0.0;
1343 flowSpeedNormal=porosity*(
n[0]*velocity_star[0] +
1344 n[1]*velocity_star[1]);
1345 if (isDOFBoundary_u != 1)
1347 dflux_mass_du +=
n[0]*df_mass_du[0];
1348 if (flowSpeedNormal < 0.0)
1349 dflux_umom_du -= flowSpeedNormal;
1353 dflux_mass_du +=
n[0]*df_mass_du[0];
1354 if (flowSpeedNormal < 0.0)
1355 dflux_umom_du -= flowSpeedNormal;
1357 if (isDOFBoundary_v != 1)
1359 dflux_mass_dv +=
n[1]*df_mass_dv[1];
1360 if (flowSpeedNormal < 0.0)
1361 dflux_vmom_dv -= flowSpeedNormal;
1365 dflux_mass_dv +=
n[1]*df_mass_dv[1];
1366 if (flowSpeedNormal < 0.0)
1367 dflux_vmom_dv -= flowSpeedNormal;
1401 if (isFluxBoundary_u == 1)
1403 dflux_umom_dp = 0.0;
1404 dflux_umom_du = 0.0;
1405 dflux_umom_dv = 0.0;
1408 if (isFluxBoundary_v == 1)
1410 dflux_vmom_dp = 0.0;
1411 dflux_vmom_du = 0.0;
1412 dflux_vmom_dv = 0.0;
1429 const int& isDOFBoundary,
1430 const int& isFluxBoundary,
1431 const double n[nSpace],
1434 const double& bc_flux,
1436 const double grad_potential[nSpace],
1438 const double& penalty,
1441 double diffusiveVelocityComponent_I,penaltyFlux,max_a;
1442 if(isFluxBoundary == 1)
1446 else if(isDOFBoundary == 1)
1450 for(
int I=0;I<nSpace;I++)
1452 diffusiveVelocityComponent_I=0.0;
1453 for(
int m=rowptr[I];m<rowptr[I+1];m++)
1455 diffusiveVelocityComponent_I -= a[m]*grad_potential[colind[m]];
1456 max_a = fmax(max_a,a[m]);
1458 flux+= diffusiveVelocityComponent_I*
n[I];
1460 penaltyFlux = max_a*penalty*(
u-bc_u);
1461 flux += penaltyFlux;
1467 std::cerr<<
"RANS3PF2D: warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl;
1477 const int& isDOFBoundary,
1478 const int& isFluxBoundary,
1479 const double n[nSpace],
1482 const double grad_v[nSpace],
1483 const double& penalty)
1485 double dvel_I,tmp=0.0,max_a=0.0;
1486 if(isFluxBoundary==0 && isDOFBoundary==1)
1488 for(
int I=0;I<nSpace;I++)
1491 for(
int m=rowptr[I];m<rowptr[I+1];m++)
1493 dvel_I -= a[m]*grad_v[colind[m]];
1494 max_a = fmax(max_a,a[m]);
1498 tmp +=max_a*penalty*
v;
1509 res[0] = grad_u[0]*
n[0]+grad_u[1]*
n[1];
1510 res[1] = grad_v[0]*
n[0]+grad_v[1]*
n[1];
1514 return u[0]*
v[1]-
u[1]*
v[0];
1518 return u[0]*
v[0]+
u[1]*
v[1];
1520 int get_distance_to_ball(
int n_balls,
double* ball_center,
double* ball_radius,
double x,
double y,
double z,
double& distance)
1525 for (
int i=0; i<n_balls; ++i)
1527 d_ball_i = std::sqrt((ball_center[i*3+0]-x)*(ball_center[i*3+0]-x)
1528 +(ball_center[i*3+1]-y)*(ball_center[i*3+1]-y)
1531 if(d_ball_i<distance)
1533 distance = d_ball_i;
1541 double x,
double y,
double z,
1544 distance = std::sqrt((ball_center[I*3+0]-x)*(ball_center[I*3+0]-x)
1545 + (ball_center[I*3+1]-y)*(ball_center[I*3+1]-y)
1551 double x,
double y,
double z,
1552 double& nx,
double& ny)
1554 double distance = std::sqrt((ball_center[I*3+0]-x)*(ball_center[I*3+0]-x)
1555 + (ball_center[I*3+1]-y)*(ball_center[I*3+1]-y)
1558 nx = (x - ball_center[I*3+0])/(distance+1e-10);
1559 ny = (y - ball_center[I*3+1])/(distance+1e-10);
1562 double* ball_velocity,
double* ball_angular_velocity,
1564 double x,
double y,
double z,
1565 double&
vx,
double&
vy)
1567 vx = ball_velocity[3*I + 0] - ball_angular_velocity[3*I + 2]*(y-ball_center[3*I + 1]);
1568 vy = ball_velocity[3*I + 1] + ball_angular_velocity[3*I + 2]*(x-ball_center[3*I + 0]);
1574 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1575 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1576 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1577 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1578 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1579 double PSTAB = args.
scalar<
double>(
"PSTAB");
1580 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1581 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
1582 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1583 int nDOF_per_element_pressure = args.
scalar<
int>(
"nDOF_per_element_pressure");
1584 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1585 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1586 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1587 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1588 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
1589 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
1590 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
1591 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
1592 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1593 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1594 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
1595 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1596 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1597 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1598 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1599 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1600 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1601 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1602 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1603 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1604 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1605 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1606 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1607 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1608 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1609 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1610 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1611 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1612 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1613 double hFactor = args.
scalar<
double>(
"hFactor");
1614 int nElements_global = args.
scalar<
int>(
"nElements_global");
1615 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1616 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
1617 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1618 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
1619 double useRBLES = args.
scalar<
double>(
"useRBLES");
1620 double useMetrics = args.
scalar<
double>(
"useMetrics");
1621 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1622 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1623 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1624 double sigma = args.
scalar<
double>(
"sigma");
1629 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1630 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1631 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1632 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1633 double C_dc = args.
scalar<
double>(
"C_dc");
1634 double C_b = args.
scalar<
double>(
"C_b");
1635 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1636 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
1637 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
1638 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
1639 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
1640 xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1641 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
1642 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
1643 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1644 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
1645 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
1646 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1647 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1648 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1649 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1650 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1651 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1652 xt::pyarray<double>& q_eddy_viscosity = args.
array<
double>(
"q_eddy_viscosity");
1653 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1654 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1655 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1656 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1657 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1658 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1659 xt::pyarray<double>& u_dof_old = args.
array<
double>(
"u_dof_old");
1660 xt::pyarray<double>& v_dof_old = args.
array<
double>(
"v_dof_old");
1661 xt::pyarray<double>& w_dof_old = args.
array<
double>(
"w_dof_old");
1662 xt::pyarray<double>& u_dof_old_old = args.
array<
double>(
"u_dof_old_old");
1663 xt::pyarray<double>& v_dof_old_old = args.
array<
double>(
"v_dof_old_old");
1664 xt::pyarray<double>& w_dof_old_old = args.
array<
double>(
"w_dof_old_old");
1665 xt::pyarray<double>& uStar_dof = args.
array<
double>(
"uStar_dof");
1666 xt::pyarray<double>& vStar_dof = args.
array<
double>(
"vStar_dof");
1667 xt::pyarray<double>& wStar_dof = args.
array<
double>(
"wStar_dof");
1668 xt::pyarray<double>& g = args.
array<
double>(
"g");
1669 const double useVF = args.
scalar<
double>(
"useVF");
1670 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1671 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1672 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
1673 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1674 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1675 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1676 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1677 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1678 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1679 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1680 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1681 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1682 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
1683 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
1684 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1685 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
1686 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1687 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1688 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1689 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1690 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1691 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1692 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1693 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1694 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1695 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1696 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1697 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1698 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1699 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1700 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1701 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1702 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1703 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1704 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1705 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1706 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1707 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1708 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1709 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1710 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1711 int offset_p = args.
scalar<
int>(
"offset_p");
1712 int offset_u = args.
scalar<
int>(
"offset_u");
1713 int offset_v = args.
scalar<
int>(
"offset_v");
1714 int offset_w = args.
scalar<
int>(
"offset_w");
1715 int stride_p = args.
scalar<
int>(
"stride_p");
1716 int stride_u = args.
scalar<
int>(
"stride_u");
1717 int stride_v = args.
scalar<
int>(
"stride_v");
1718 int stride_w = args.
scalar<
int>(
"stride_w");
1719 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1720 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1721 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1722 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
1723 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1724 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1725 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1726 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1727 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1728 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1729 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1730 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1731 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
1732 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1733 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1734 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1735 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1736 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1737 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1738 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1739 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1740 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1741 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1742 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1743 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1744 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1745 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1746 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1747 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1748 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1749 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1750 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1751 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1752 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1753 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1754 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1755 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1756 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1757 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1758 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1759 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1760 xt::pyarray<double>& q_grad_u = args.
array<
double>(
"q_grad_u");
1761 xt::pyarray<double>& q_grad_v = args.
array<
double>(
"q_grad_v");
1762 xt::pyarray<double>& q_grad_w = args.
array<
double>(
"q_grad_w");
1763 xt::pyarray<double>& q_divU = args.
array<
double>(
"q_divU");
1764 xt::pyarray<double>& ebqe_grad_u = args.
array<
double>(
"ebqe_grad_u");
1765 xt::pyarray<double>& ebqe_grad_v = args.
array<
double>(
"ebqe_grad_v");
1766 xt::pyarray<double>& ebqe_grad_w = args.
array<
double>(
"ebqe_grad_w");
1767 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1768 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1769 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
1770 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1771 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1772 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1773 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1774 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1775 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1776 xt::pyarray<double>& q_rho = args.
array<
double>(
"q_rho");
1777 xt::pyarray<double>& ebqe_rho = args.
array<
double>(
"ebqe_rho");
1778 xt::pyarray<double>& q_nu = args.
array<
double>(
"q_nu");
1779 xt::pyarray<double>& ebqe_nu = args.
array<
double>(
"ebqe_nu");
1780 int nParticles = args.
scalar<
int>(
"nParticles");
1781 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
1782 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
1783 double particle_beta = args.
scalar<
double>(
"particle_beta");
1784 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
1785 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
1786 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
1787 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
1788 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
1789 xt::pyarray<double>& particle_netForces = args.
array<
double>(
"particle_netForces");
1790 xt::pyarray<double>& particle_netMoments = args.
array<
double>(
"particle_netMoments");
1791 xt::pyarray<double>& particle_surfaceArea = args.
array<
double>(
"particle_surfaceArea");
1792 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
1793 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
1794 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
1795 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
1796 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
1797 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
1798 xt::pyarray<double>& phisError = args.
array<
double>(
"phisError");
1799 xt::pyarray<double>& phisErrorNodal = args.
array<
double>(
"phisErrorNodal");
1800 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
1801 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
1803 double cE = args.
scalar<
double>(
"cE");
1804 int MULTIPLY_EXTERNAL_FORCE_BY_DENSITY = args.
scalar<
int>(
"MULTIPLY_EXTERNAL_FORCE_BY_DENSITY");
1805 xt::pyarray<double>& forcex = args.
array<
double>(
"forcex");
1806 xt::pyarray<double>& forcey = args.
array<
double>(
"forcey");
1807 xt::pyarray<double>& forcez = args.
array<
double>(
"forcez");
1808 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
1809 double dt = args.
scalar<
double>(
"dt");
1810 xt::pyarray<double>& quantDOFs = args.
array<
double>(
"quantDOFs");
1811 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
1812 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
1813 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
1814 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
1815 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
1816 double order_polynomial = args.
scalar<
double>(
"order_polynomial");
1817 xt::pyarray<double>& isActiveDOF = args.
array<
double>(
"isActiveDOF");
1818 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
1819 xt::pyarray<double>& ncDrag = args.
array<
double>(
"ncDrag");
1820 xt::pyarray<double>& betaDrag = args.
array<
double>(
"betaDrag");
1821 xt::pyarray<double>& vos_vel_nodes = args.
array<
double>(
"vos_vel_nodes");
1822 xt::pyarray<double>& entropyResidualPerNode = args.
array<
double>(
"entropyResidualPerNode");
1823 xt::pyarray<double>& laggedEntropyResidualPerNode = args.
array<
double>(
"laggedEntropyResidualPerNode");
1824 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
1825 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
1826 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
1827 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
1828 int NNZ_1D = args.
scalar<
int>(
"NNZ_1D");
1829 xt::pyarray<int>& csrRowIndeces_1D = args.
array<
int>(
"csrRowIndeces_1D");
1830 xt::pyarray<int>& csrColumnOffsets_1D = args.
array<
int>(
"csrColumnOffsets_1D");
1831 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
1832 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
1833 xt::pyarray<double>& isBoundary_1D = args.
array<
double>(
"isBoundary_1D");
1834 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
1835 gf.useExact=useExact;
1836 gf_s.useExact=useExact;
1840 double cut_cell_boundary_length=0.0, p_force_x=0.0, p_force_y=0.0;
1841 register double element_uStar_He[nElements_global], element_vStar_He[nElements_global];
1844 den_hi.resize(numDOFs_1D,0.0);
1854 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1860 if (
psi.size() != numDOFs_1D)
1861 psi.resize(numDOFs_1D);
1862 for (
int i=0; i<NNZ_1D; i++)
1864 uStar_dMatrix[i]=0.;
1865 vStar_dMatrix[i]=0.;
1869 for (
int i=0; i<numDOFs_1D; i++)
1873 entropyResidualPerNode[i]=0.;
1883 double mesh_volume_conservation=0.0,
1884 mesh_volume_conservation_weak=0.0,
1885 mesh_volume_conservation_err_max=0.0,
1886 mesh_volume_conservation_err_max_weak=0.0;
1887 double globalConservationError=0.0;
1888 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
1889 for(
int eN=0;eN<nElements_global;eN++)
1891 register double elementTransport[nDOF_test_element][nDOF_trial_element];
1892 register double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1894 register double elementResidual_p[nDOF_test_element],elementResidual_mesh[nDOF_test_element],
1895 elementResidual_u[nDOF_test_element],
1896 elementResidual_v[nDOF_test_element],
1897 mom_u_source_i[nDOF_test_element],
1898 mom_v_source_i[nDOF_test_element],
1899 betaDrag_i[nDOF_test_element],
1900 vos_i[nDOF_test_element],
1901 phisErrorElement[nDOF_test_element],
1903 elementEntropyResidual[nDOF_test_element],
1906 double element_active=1.0;
1907 double mesh_volume_conservation_element=0.0,
1908 mesh_volume_conservation_element_weak=0.0;
1910 double linVisc_eN = 0, nlinVisc_eN_num = 0, nlinVisc_eN_den = 0;
1912 double det_hess_uStar_Ke=0.0, det_hess_vStar_Ke=0.0, area_Ke=0.0;
1913 for (
int i=0;i<nDOF_test_element;i++)
1915 int eN_i = eN*nDOF_test_element+i;
1916 elementResidual_p_save[eN_i]=0.0;
1917 elementResidual_mesh[i]=0.0;
1918 elementResidual_p[i]=0.0;
1919 elementResidual_u[i]=0.0;
1920 elementResidual_v[i]=0.0;
1921 mom_u_source_i[i]=0.0;
1922 mom_v_source_i[i]=0.0;
1925 phisErrorElement[i]=0.0;
1927 elementEntropyResidual[i]=0.0;
1928 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1930 for (
int j=0;j<nDOF_trial_element;j++)
1932 elementTransport[i][j]=0.0;
1933 elementTransposeTransport[i][j]=0.0;
1938 if(use_ball_as_particle==1)
1940 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1942 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1943 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1944 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1945 phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]]);
1952 double _distance[nDOF_mesh_trial_element]={0.0};
1954 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1956 if(use_ball_as_particle==1)
1959 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1960 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1961 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1966 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
1968 if ( _distance[I] > 0)
1971 if (pos_counter == 3)
1973 element_active = 1.0;
1975 else if (pos_counter == 0)
1977 element_active = 1.0;
1981 element_active = 1.0;
1984 double sub_mesh_dof[6*3], sub_u_dof[15], sub_v_dof[15], sub_phi_dof[6], sub_p_dof[6];
1985 int boundaryNodes[6] = {0,0,0,0,0,0};
1986 std::vector<int> ls_nodes;
1987 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1989 for (
int K=0;K<nDOF_mesh_trial_element;K++)
1997 const double eps = 1.0e-4;
1998 double delta_phi=0.0,theta;
1999 delta_phi = _distance[(I+1)%3] - _distance[I];
2000 if (fabs(delta_phi) > eps)
2003 theta = -_distance[I]/delta_phi;
2004 if (theta > 1.0-eps || theta < eps)
2006 if (theta > 1.0-eps && theta <= 1.0)
2008 ls_nodes.push_back((I+1)%3);
2012 else if (theta > 0.0 && theta < eps)
2014 ls_nodes.push_back(I);
2022 boundaryNodes[3+I]=1;
2023 ls_nodes.push_back(3+I);
2029 if (fabs(_distance[I]) <= eps)
2032 boundaryNodes[3+I]=1;
2033 boundaryNodes[(I+1)%3]=1;
2034 ls_nodes.push_back(I);
2035 ls_nodes.push_back((I+1)%3);
2038 assert(theta <= 1.0);
2039 GI[3*3 + I*3 + I] = 1.0-theta;
2040 GI[3*3 + I*3 + (I+1)%3] = theta;
2041 GI[3*3 + I*3 + (I+2)%3] = 0.0;
2043 if (ls_nodes.size() != 2)
2045 std::cout<<
"level set nodes not 2 "<<ls_nodes.size()<<std::endl;
2046 for(
int i=0;i<ls_nodes.size();i++)
2047 std::cout<<ls_nodes[i]<<std::endl;
2048 std::sort(ls_nodes.begin(),ls_nodes.end());
2050 int sub_mesh_l2g[12] = {0,3,5,
2054 for (
int I=0; I<6; I++)
2056 sub_phi_dof[I] = 0.0;
2058 for (
int K=0; K<3; K++)
2059 sub_mesh_dof[I*3+K] = 0.0;
2060 for (
int J=0; J<3; J++)
2062 for (
int K=0; K<3; K++)
2064 sub_mesh_dof[I*3+K] += GI[I*3+J]*mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+J]+K];
2066 sub_phi_dof[I] += GI[I*3+J]*phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+J]];
2067 sub_p_dof[I] += GI[I*3+J]*p_dof[p_l2g[eN*nDOF_per_element_pressure+J]];
2070 int L = ls_nodes[0],
R=ls_nodes[1];
2071 double DX=sub_mesh_dof[
L*3+0] - sub_mesh_dof[
R*3+0];
2072 double DY=sub_mesh_dof[
L*3+1] - sub_mesh_dof[
R*3+1];
2073 double DS = std::sqrt(DX*DX+DY*DY);
2074 double nx = -DY/DS, ny = DX/DS;
2075 double nxL,nyL,nxR,nyR;
2078 sub_mesh_dof[
L*3+0],sub_mesh_dof[
L*3+1],0.0,
2082 sub_mesh_dof[
R*3+0],sub_mesh_dof[
R*3+1],0.0,
2086 double n_fluid_sign = -
sgn(nx*0.5*(nxL+nxR)+ny*0.5*(nyL+nyR));
2091 cut_cell_boundary_length += DS;
2092 p_force_x += sub_p_dof[
L]*nx*0.5*DS + sub_p_dof[
R]*nx*0.5*DS;
2093 p_force_y += sub_p_dof[
L]*ny*0.5*DS + sub_p_dof[
R]*ny*0.5*DS;
2100 double lagrangeNodes[9*3];
2101 for (
int K=0;K<3;K++)
2103 lagrangeNodes[0*3+K] = 0.5*(sub_mesh_dof[0*3+K] + sub_mesh_dof[3*3+0*3+K]);
2104 lagrangeNodes[1*3+K] = 0.5*(sub_mesh_dof[1*3+K] + sub_mesh_dof[3*3+0*3+K]);
2105 lagrangeNodes[2*3+K] = 0.5*(sub_mesh_dof[1*3+K] + sub_mesh_dof[3*3+1*3+K]);
2106 lagrangeNodes[3*3+K] = 0.5*(sub_mesh_dof[2*3+K] + sub_mesh_dof[3*3+1*3+K]);
2107 lagrangeNodes[4*3+K] = 0.5*(sub_mesh_dof[2*3+K] + sub_mesh_dof[3*3+2*3+K]);
2108 lagrangeNodes[5*3+K] = 0.5*(sub_mesh_dof[0*3+K] + sub_mesh_dof[3*3+2*3+K]);
2109 lagrangeNodes[6*3+K] = 0.5*(sub_mesh_dof[3*3+0*3+K] + sub_mesh_dof[3*3+1*3+K]);
2110 lagrangeNodes[7*3+K] = 0.5*(sub_mesh_dof[3*3+1*3+K] + sub_mesh_dof[3*3+2*3+K]);
2111 lagrangeNodes[8*3+K] = 0.5*(sub_mesh_dof[3*3+2*3+K] + sub_mesh_dof[3*3+0*3+K]);
2114 for (
int I=0;I<6;I++)
2116 baryCoords(&sub_mesh_dof[0],&sub_mesh_dof[1*3],&sub_mesh_dof[2*3],&sub_mesh_dof[I*3],lambda);
2118 G2I[I*6+0] = lambda[0]*(2.0*lambda[0] - 1.0);
2119 G2I[I*6+1] = lambda[1]*(2.0*lambda[1] - 1.0);
2120 G2I[I*6+2] = lambda[2]*(2.0*lambda[2] - 1.0);
2121 G2I[I*6+3] = 4.0*lambda[0]*lambda[1];
2122 G2I[I*6+4] = 4.0*lambda[1]*lambda[2];
2123 G2I[I*6+5] = 4.0*lambda[2]*lambda[0];
2125 for (
int I=0;I<9;I++)
2127 baryCoords(&sub_mesh_dof[0],&sub_mesh_dof[1*3],&sub_mesh_dof[2*3],&lagrangeNodes[I*3],lambda);
2128 G2I[6*6 + I*6 + 0] = lambda[0]*(2.0*lambda[0] - 1.0);
2129 G2I[6*6 + I*6 + 1] = lambda[1]*(2.0*lambda[1] - 1.0);
2130 G2I[6*6 + I*6 + 2] = lambda[2]*(2.0*lambda[2] - 1.0);
2131 G2I[6*6 + I*6 + 3] = 4.0*lambda[0]*lambda[1];
2132 G2I[6*6 + I*6 + 4] = 4.0*lambda[1]*lambda[2];
2133 G2I[6*6 + I*6 + 5] = 4.0*lambda[2]*lambda[0];
2135 for (
int I=0; I<15; I++)
2139 for (
int J=0; J<6; J++)
2141 sub_u_dof[I] += G2I[I*6+J]*u_dof[vel_l2g[eN*nDOF_trial_element+J]];
2142 sub_v_dof[I] += G2I[I*6+J]*v_dof[vel_l2g[eN*nDOF_trial_element+J]];
2145 for (
int esN=0;esN<4;esN++)
2147 std::cout<<sub_mesh_l2g[esN*3]<<
'\t'<<sub_mesh_l2g[esN*3+1]<<
'\t'<<sub_mesh_l2g[esN*3+2]<<std::endl;
2149 for (
int I=0; I<6; I++)
2151 std::cout<<sub_mesh_dof[I*3+0]<<
'\t'<<sub_mesh_dof[I*3+1]<<
'\t'<<sub_mesh_dof[I*3+2]<<
'\t'<<boundaryNodes[I]<<
'\t'<<sub_phi_dof[I]<<
'\t'<<sub_p_dof[I]<<
'\t'<<sub_u_dof[I]<<
'\t'<<sub_v_dof[I]<<
'\t'<<G2I[I*6+0]<<
'\t'<<G2I[I*6+1]<<
'\t'<<G2I[I*6+2]<<
'\t'<<G2I[I*6+3]<<
'\t'<<G2I[I*6+4]<<
'\t'<<G2I[I*6+5]<<std::endl;
2160 double _distance[nDOF_mesh_trial_element]={0.0};
2162 for (
int I=0;I<nDOF_mesh_trial_element;I++)
2164 if(use_ball_as_particle==1)
2167 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
2168 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
2169 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
2174 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
2176 if ( _distance[I] >= 0)
2179 if (pos_counter == 2)
2183 for (
int I=0;I<nDOF_mesh_trial_element;I++)
2186 if (_distance[I] < 0)
2192 assert(opp_node >=0);
2193 assert(opp_node <nDOF_mesh_trial_element);
2198 const int ebN = elementBoundariesArray[eN*nDOF_mesh_trial_element+opp_node];
2199 const int eN_oppo = (eN == elementBoundaryElementsArray[ebN*2+0])?elementBoundaryElementsArray[ebN*2+1]:elementBoundaryElementsArray[ebN*2+0];
2200 if((mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%3]<nNodes_owned
2201 || mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%3]<nNodes_owned)
2207 if (eN == elementBoundaryElementsArray[ebN*2+0])
2214 if(use_ball_as_particle==1)
2216 double middle_point_coord[3]={0.0};
2217 double middle_point_distance;
2218 middle_point_coord[0] = 0.5*(mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%3]+0]+mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%3]+0]);
2219 middle_point_coord[1] = 0.5*(mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%3]+1]+mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%3]+1]);
2221 middle_point_coord[0],middle_point_coord[1],middle_point_coord[2],
2222 middle_point_distance);
2232 double distance=1e10, distance_to_ith_particle;
2233 for (
int i=0;i<nParticles;++i)
2235 distance_to_ith_particle=particle_signed_distances[i*nElements_global*nQuadraturePoints_element
2236 +eN*nQuadraturePoints_element
2238 if (distance_to_ith_particle<distance)
2240 distance = distance_to_ith_particle;
2250 if(ebN<nElementBoundaries_owned)
2252 assert(eN_oppo==-1);
2256 else if (pos_counter == 3)
2259 for (
int i=0;i<nDOF_test_element;i++)
2261 isActiveDOF[offset_u+stride_u*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
2262 isActiveDOF[offset_v+stride_v*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
2270 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
2271 for (
int j=0;j<nDOF_mesh_trial_element;j++)
2273 register int eN_j = eN*nDOF_mesh_trial_element+j;
2274 element_phi[j] = phi_dof[p_l2g[eN_j]];
2275 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
2277 double element_nodes[nDOF_mesh_trial_element*3];
2278 for (
int i=0;i<nDOF_mesh_trial_element;i++)
2280 register int eN_i=eN*nDOF_mesh_trial_element+i;
2281 for(
int I=0;I<3;I++)
2282 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
2284 gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),
false);
2285 gf.calculate(element_phi, element_nodes, x_ref.data(),
false);
2289 for(
int k=0;k<nQuadraturePoints_element;k++)
2294 register int eN_k = eN*nQuadraturePoints_element+k,
2295 eN_k_nSpace = eN_k*nSpace,
2297 eN_nDOF_trial_element = eN*nDOF_trial_element;
2298 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,un=0.0,vn=0.0,wn=0.0,
2299 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
2308 dmass_adv_u[nSpace],
2309 dmass_adv_v[nSpace],
2310 dmass_adv_w[nSpace],
2312 dmom_u_adv_u[nSpace],
2313 dmom_u_adv_v[nSpace],
2314 dmom_u_adv_w[nSpace],
2316 dmom_v_adv_u[nSpace],
2317 dmom_v_adv_v[nSpace],
2318 dmom_v_adv_w[nSpace],
2320 dmom_w_adv_u[nSpace],
2321 dmom_w_adv_v[nSpace],
2322 dmom_w_adv_w[nSpace],
2323 mom_uu_diff_ten[nSpace],
2324 mom_vv_diff_ten[nSpace],
2325 mom_ww_diff_ten[nSpace],
2336 dmom_u_ham_grad_p[nSpace],
2337 dmom_u_ham_grad_u[nSpace],
2339 dmom_v_ham_grad_p[nSpace],
2340 dmom_v_ham_grad_v[nSpace],
2342 dmom_w_ham_grad_p[nSpace],
2343 dmom_w_ham_grad_w[nSpace],
2354 Lstar_u_p[nDOF_test_element],
2355 Lstar_v_p[nDOF_test_element],
2356 Lstar_w_p[nDOF_test_element],
2357 Lstar_u_u[nDOF_test_element],
2358 Lstar_v_v[nDOF_test_element],
2359 Lstar_w_w[nDOF_test_element],
2360 Lstar_p_u[nDOF_test_element],
2361 Lstar_p_v[nDOF_test_element],
2362 Lstar_p_w[nDOF_test_element],
2367 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
2368 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2371 jacInv[nSpace*nSpace],
2372 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
2373 vel_hess_trial[nDOF_trial_element*
nSpace2],
2374 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
2375 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
2376 u_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2377 v_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2383 dmom_u_source[nSpace],
2384 dmom_v_source[nSpace],
2385 dmom_w_source[nSpace],
2389 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
2391 ck.calculateMapping_element(eN,
2395 mesh_trial_ref.data(),
2396 mesh_grad_trial_ref.data(),
2401 ck.calculateH_element(eN,
2403 nodeDiametersArray.data(),
2405 mesh_trial_ref.data(),
2407 ck.calculateMappingVelocity_element(eN,
2409 mesh_velocity_dof.data(),
2411 mesh_trial_ref.data(),
2416 dV = fabs(jacDet)*dV_ref[k];
2417 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2420 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2421 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2422 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2426 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
2427 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
2432 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
2433 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
2436 ck.valFromDOF(u_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],un);
2437 ck.valFromDOF(v_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],vn);
2441 for (
int I=0;I<nSpace;I++)
2442 grad_p[I] = q_grad_p[eN_k_nSpace + I];
2443 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
2444 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
2445 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
2446 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
2447 ck.hessFromDOF(uStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_uStar);
2448 ck.hessFromDOF(vStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_vStar);
2451 for (
int j=0;j<nDOF_trial_element;j++)
2454 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
2455 for (
int I=0;I<nSpace;I++)
2458 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
2459 if (ARTIFICIAL_VISCOSITY==4)
2462 u_times_vel_grad_test_dV[j*nSpace+I] =
2463 u*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_u[I];
2464 v_times_vel_grad_test_dV[j*nSpace+I] =
2465 v*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_v[I];
2472 if (ARTIFICIAL_VISCOSITY==3)
2474 det_hess_uStar_Ke += (hess_uStar[0]*hess_uStar[3] - hess_uStar[2]*hess_uStar[1])*dV;
2475 det_hess_vStar_Ke += (hess_vStar[0]*hess_vStar[3] - hess_vStar[2]*hess_vStar[1])*dV;
2479 double div_mesh_velocity=0.0;
2480 int NDOF_MESH_TRIAL_ELEMENT=3;
2481 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
2483 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
2484 div_mesh_velocity +=
2485 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*nSpace+0] +
2486 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*nSpace+1];
2488 mesh_volume_conservation_element += (alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity)*dV;
2489 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
2491 porosity = 1.0 - q_vos[eN_k];
2500 double distance_to_omega_solid = 1e10;
2501 if(use_ball_as_particle==1)
2505 distance_to_omega_solid);
2509 for (
int i = 0; i < nParticles; i++)
2511 double distance_to_i_th_solid = particle_signed_distances[i * nElements_global * nQuadraturePoints_element + eN_k];
2512 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
2515 phi_solid[eN_k] = distance_to_omega_solid;
2527 elementDiameter[eN],
2528 smagorinskyConstant,
2529 turbulenceClosureModel,
2534 &normal_phi[eN_k_nSpace],
2535 distance_to_omega_solid,
2548 q_velocity_sge[eN_k_nSpace+0],
2549 q_velocity_sge[eN_k_nSpace+1],
2550 q_velocity_sge[eN_k_nSpace+1],
2551 q_eddy_viscosity[eN_k],
2598 MULTIPLY_EXTERNAL_FORCE_BY_DENSITY,
2602 MATERIAL_PARAMETERS_AS_FUNCTION,
2603 density_as_function[eN_k],
2604 dynamic_viscosity_as_function[eN_k],
2607 use_ball_as_particle,
2610 ball_velocity.data(),
2611 ball_angular_velocity.data(),
2612 INT_BY_PARTS_PRESSURE);
2615 mass_source = q_mass_source[eN_k];
2616 for (
int I=0;I<nSpace;I++)
2618 dmom_u_source[I] = 0.0;
2619 dmom_v_source[I] = 0.0;
2620 dmom_w_source[I] = 0.0;
2631 q_eddy_viscosity[eN_k],
2638 q_velocity_sge[eN_k_nSpace+0],
2639 q_velocity_sge[eN_k_nSpace+1],
2640 q_velocity_sge[eN_k_nSpace+1],
2641 eps_solid[elementFlags[eN]],
2643 q_velocity_solid[eN_k_nSpace+0],
2644 q_velocity_solid[eN_k_nSpace+1],
2645 q_velocity_solid[eN_k_nSpace+1],
2646 q_velocityStar_solid[eN_k_nSpace+0],
2647 q_velocityStar_solid[eN_k_nSpace+1],
2648 q_velocityStar_solid[eN_k_nSpace+1],
2655 q_grad_vos[eN_k_nSpace+0],
2656 q_grad_vos[eN_k_nSpace+1],
2657 q_grad_vos[eN_k_nSpace+1]);
2658 double C_particles=0.0;
2659 if(nParticles > 0 && USE_SBM==0)
2664 nQuadraturePoints_global,
2665 &particle_signed_distances[eN_k],
2666 &particle_signed_distance_normals[eN_k_3d],
2667 &particle_velocities[eN_k_3d],
2668 particle_centroids.data(),
2669 use_ball_as_particle,
2672 ball_velocity.data(),
2673 ball_angular_velocity.data(),
2675 particle_penalty_constant/h_phi,
2676 particle_alpha/h_phi,
2677 particle_beta/h_phi,
2694 q_velocity_sge[eN_k_nSpace+0],
2695 q_velocity_sge[eN_k_nSpace+1],
2696 q_velocity_sge[eN_k_nSpace+1],
2719 particle_netForces.data(),
2720 particle_netMoments.data(),
2721 particle_surfaceArea.data());
2726 nQuadraturePoints_global,
2727 &particle_signed_distances[eN_k],
2728 &particle_signed_distance_normals[eN_k_3d],
2729 &particle_velocities[eN_k_3d],
2730 particle_centroids.data(),
2731 use_ball_as_particle,
2734 ball_velocity.data(),
2735 ball_angular_velocity.data(),
2736 particle_penalty_constant/h_phi,
2737 particle_alpha/h_phi,
2738 particle_beta/h_phi,
2755 q_velocity_sge[eN_k_nSpace+0],
2756 q_velocity_sge[eN_k_nSpace+1],
2757 q_velocity_sge[eN_k_nSpace+1],
2762 particle_netForces.data(),
2763 particle_netMoments.data());
2765 if (turbulenceClosureModel >= 3)
2767 const double c_mu = 0.09;
2782 &q_turb_var_grad_0[eN_k_nSpace],
2783 q_eddy_viscosity[eN_k],
2801 q_mom_u_acc[eN_k] = mom_u_acc;
2802 q_mom_v_acc[eN_k] = mom_v_acc;
2805 q_mass_adv[eN_k_nSpace+0] =
u;
2806 q_mass_adv[eN_k_nSpace+1] =
v;
2811 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
2812 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
2814 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
2815 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
2818 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
2819 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
2821 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
2822 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
2834 if (q_dV_last[eN_k] <= -100)
2835 q_dV_last[eN_k] = dV;
2838 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2844 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2858 mom_u_acc_t *= dmom_u_acc_u;
2859 mom_v_acc_t *= dmom_v_acc_v;
2865 ck.Mass_strong(-q_dvos_dt[eN_k]) +
2866 ck.Advection_strong(dmass_adv_u,grad_u) +
2867 ck.Advection_strong(dmass_adv_v,grad_v) +
2869 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
2871 ck.Reaction_strong(mass_source);
2874 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
2875 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
2879 ck.Mass_strong(mom_u_acc_t) +
2880 ck.Advection_strong(dmom_adv_sge,grad_u) +
2881 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
2882 ck.Reaction_strong(mom_u_source) -
2883 ck.Reaction_strong(
u*div_mesh_velocity);
2886 ck.Mass_strong(mom_v_acc_t) +
2887 ck.Advection_strong(dmom_adv_sge,grad_v) +
2888 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
2889 ck.Reaction_strong(mom_v_source) -
2890 ck.Reaction_strong(
v*div_mesh_velocity);
2900 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
2902 elementDiameter[eN],
2907 dmom_u_ham_grad_p[0],
2917 dmom_u_ham_grad_p[0],
2922 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
2923 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
2936 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
2937 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
2940 mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge[eN_k_nSpace+0]);
2941 mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge[eN_k_nSpace+0]);
2945 for (
int i=0;i<nDOF_test_element;i++)
2947 register int i_nSpace = i*nSpace;
2952 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2953 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2955 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2956 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2960 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2961 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2966 if (ARTIFICIAL_VISCOSITY==0 || ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2968 q_numDiff_u[eN_k] = 0;
2969 q_numDiff_v[eN_k] = 0;
2970 q_numDiff_w[eN_k] = 0;
2972 else if (ARTIFICIAL_VISCOSITY==1)
2974 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v);
2975 q_numDiff_u[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2976 (1.0-useMetrics)*hFactor*hFactor*elementDiameter[eN]*elementDiameter[eN]);
2977 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
2978 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
2982 double rho = q_rho[eN_k];
2983 double mu = q_rho[eN_k]*q_nu[eN_k];
2985 double vel2 =
u*
u +
v*
v;
2989 porosity*rho*((
u-un)/dt + (
u*grad_u[0]+
v*grad_u[1]) - g[0])
2990 + (KILL_PRESSURE_TERM == 1 ? 0. : 1.)*grad_p[0]
2991 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcex[eN_k]
2992 - mu*(hess_u[0] + hess_u[3])
2993 - mu*(hess_u[0] + hess_v[2]);
2996 porosity*rho*((
v-vn)/dt + (
u*grad_v[0]+
v*grad_v[1]) - g[1])
2997 + (KILL_PRESSURE_TERM == 1 ? 0. : 1.)*grad_p[1]
2998 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcey[eN_k]
2999 - mu*(hess_v[0] + hess_v[3])
3000 - mu*(hess_u[1] + hess_v[3]);
3003 double entRes_times_u = Res_in_x*
u + Res_in_y*
v;
3005 double hK = elementDiameter[eN]/order_polynomial;
3006 q_numDiff_u[eN_k] = fmin(
cMax*porosity*rho*hK*std::sqrt(vel2),
3007 cE*hK*hK*fabs(entRes_times_u)/(vel2+1E-10));
3008 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
3009 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
3013 linVisc_eN = fmax(porosity*rho*std::sqrt(vel2),linVisc_eN);
3014 nlinVisc_eN_num = fmax(fabs(entRes_times_u),nlinVisc_eN_num);
3015 nlinVisc_eN_den = fmax(vel2,nlinVisc_eN_den);
3026 q_velocity[eN_k_nSpace+0]=
u;
3027 q_velocity[eN_k_nSpace+1]=
v;
3029 for (
int I=0;I<nSpace;I++)
3031 q_grad_u[eN_k_nSpace+I] = grad_u[I];
3032 q_grad_v[eN_k_nSpace+I] = grad_v[I];
3036 q_divU[eN_k] = q_grad_u[eN_k_nSpace+0] + q_grad_v[eN_k_nSpace+1];
3039 double unit_normal[nSpace];
3040 double norm_grad_phi = 0.;
3041 for (
int I=0;I<nSpace;I++)
3042 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
3043 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
3044 for (
int I=0;I<nSpace;I++)
3045 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
3049 v1[0]=1.-unit_normal[0]*unit_normal[0];
3050 v1[1]=-unit_normal[0]*unit_normal[1];
3053 v2[0]=-unit_normal[0]*unit_normal[1];
3054 v2[1]=1.-unit_normal[1]*unit_normal[1];
3055 double delta =
gf.D(eps_mu,
phi[eN_k]);
3056 register double vel_tgrad_test_i[nSpace], tgrad_u[nSpace], tgrad_v[nSpace];
3065 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3067 velStar[0] = q_velocity_sge[eN_k_nSpace+0];
3068 velStar[1] = q_velocity_sge[eN_k_nSpace+1];
3071 for(
int i=0;i<nDOF_test_element;i++)
3073 register int i_nSpace=i*nSpace;
3075 &vel_grad_trial[i_nSpace],
3077 phisErrorElement[i]+=std::abs(phisError[eN_k_nSpace+0])*p_test_dV[i];
3095 elementResidual_u[i] +=
3096 ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
3097 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
3098 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
3099 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
3101 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
3102 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
3103 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+0] : 0.) +
3105 USE_SUPG*
ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
3106 ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]) +
3108 ck.NumericalDiffusion(delta*sigma*dV,v1,vel_tgrad_test_i) +
3109 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_u,vel_tgrad_test_i);
3110 mom_u_source_i[i] +=
ck.Reaction_weak(mom_u_source,vel_test_dV[i]);
3111 betaDrag_i[i] +=
ck.Reaction_weak(dmom_u_source[0],
3113 vos_i[i] +=
ck.Reaction_weak(1.0-porosity,
3116 elementResidual_v[i] +=
3117 ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
3118 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
3119 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]) +
3120 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
3122 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
3123 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
3124 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+1] : 0.) +
3126 USE_SUPG*
ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
3127 ck.NumericalDiffusion(q_numDiff_v_last[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]) +
3129 ck.NumericalDiffusion(delta*sigma*dV,v2,vel_tgrad_test_i) +
3130 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_v,vel_tgrad_test_i);
3131 mom_v_source_i[i] +=
ck.Reaction_weak(mom_v_source,vel_test_dV[i]);
3145 if (ARTIFICIAL_VISCOSITY==4)
3149 elementEntropyResidual[i] +=
3151 ck.Mass_weak(mom_u_acc_t,
u*vel_test_dV[i]) +
3152 ck.Advection_weak(mom_u_adv,&u_times_vel_grad_test_dV[i_nSpace])+
3153 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),
3154 sdInfo_u_u_colind.data(),
3157 &u_times_vel_grad_test_dV[i_nSpace]) +
3158 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),
3159 sdInfo_u_v_colind.data(),
3162 &u_times_vel_grad_test_dV[i_nSpace]) +
3163 ck.Reaction_weak(mom_u_source,
u*vel_test_dV[i]) +
3164 ck.Hamiltonian_weak(mom_u_ham,
u*vel_test_dV[i])
3166 ck.Mass_weak(mom_v_acc_t,
v*vel_test_dV[i]) +
3167 ck.Advection_weak(mom_v_adv,&v_times_vel_grad_test_dV[i_nSpace])+
3168 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),
3169 sdInfo_v_u_colind.data(),
3172 &v_times_vel_grad_test_dV[i_nSpace])+
3173 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),
3174 sdInfo_v_v_colind.data(),
3177 &v_times_vel_grad_test_dV[i_nSpace])+
3178 ck.Reaction_weak(mom_v_source,
v*vel_test_dV[i]) +
3179 ck.Hamiltonian_weak(mom_v_ham,
v*vel_test_dV[i]);
3181 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3183 for(
int j=0;j<nDOF_trial_element;j++)
3185 int j_nSpace = j*nSpace;
3186 int i_nSpace = i*nSpace;
3187 elementTransport[i][j] +=
3188 q_rho[eN_k]*porosity*
3189 ck.AdvectionJacobian_strong(velStar,
3190 &vel_grad_test_dV[j_nSpace])
3191 *vel_trial_ref[k*nDOF_trial_element+i];
3192 elementTransposeTransport[i][j] +=
3193 q_rho[eN_k]*porosity*
3194 ck.AdvectionJacobian_strong(velStar,
3195 &vel_grad_test_dV[i_nSpace])
3196 *vel_trial_ref[k*nDOF_trial_element+j];
3201 element_uStar_He[eN] = det_hess_uStar_Ke/area_Ke;
3202 element_vStar_He[eN] = det_hess_vStar_Ke/area_Ke;
3207 double hK = elementDiameter[eN];
3208 double artVisc = fmin(
cMax*hK*linVisc_eN,
3209 cE*hK*hK*nlinVisc_eN_num/(nlinVisc_eN_den+1E-10));
3210 for(
int k=0;k<nQuadraturePoints_element;k++)
3212 register int eN_k = eN*nQuadraturePoints_element+k;
3213 q_numDiff_u[eN_k] = artVisc;
3214 q_numDiff_v[eN_k] = artVisc;
3215 q_numDiff_w[eN_k] = artVisc;
3221 for(
int i=0;i<nDOF_test_element;i++)
3223 register int eN_i=eN*nDOF_test_element+i;
3224 phisErrorNodal[vel_l2g[eN_i]]+= element_active*phisErrorElement[i];
3228 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=element_active*elementResidual_u[i];
3229 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=element_active*elementResidual_v[i];
3231 ncDrag[offset_u+stride_u*vel_l2g[eN_i]]+=mom_u_source_i[i];
3232 ncDrag[offset_v+stride_v*vel_l2g[eN_i]]+=mom_v_source_i[i];
3233 betaDrag[vel_l2g[eN_i]] += betaDrag_i[i];
3234 vos_vel_nodes[vel_l2g[eN_i]] += vos_i[i];
3237 if (ARTIFICIAL_VISCOSITY==3)
3239 uStar_hi[vel_l2g[eN_i]] += element_uStar_He[eN];
3240 vStar_hi[vel_l2g[eN_i]] += element_vStar_He[eN];
3241 den_hi[vel_l2g[eN_i]] += 1;
3243 if (ARTIFICIAL_VISCOSITY==4)
3246 entropyResidualPerNode[vel_l2g[eN_i]] += elementEntropyResidual[i];
3248 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3250 for (
int j=0;j<nDOF_trial_element;j++)
3252 int eN_i_j = eN_i*nDOF_trial_element+j;
3254 + csrColumnOffsets_1D[eN_i_j]]
3255 += elementTransport[i][j];
3258 + csrColumnOffsets_1D[eN_i_j]]
3259 += elementTransposeTransport[i][j];
3270 std::cout<<std::flush;
3272 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3275 for (
int i=0; i<numDOFs_1D; i++)
3277 if (ARTIFICIAL_VISCOSITY==4)
3280 double max_u2i = (std::pow(u_dof[i],2.) +
3281 std::pow(v_dof[i],2.));
3282 double min_u2i = max_u2i;
3283 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3285 int j = colind_1D[offset];
3286 double u2j = (std::pow(u_dof[j],2.) +
3287 std::pow(v_dof[j],2.));
3288 max_u2i = fmax(max_u2i,u2j);
3289 min_u2i = fmin(min_u2i,u2j);
3291 double normi = 0.5*(max_u2i + min_u2i) + 1E-10;
3292 entropyResidualPerNode[i] = fabs(entropyResidualPerNode[i])/normi;
3297 double uStari = uStar_dof[i];
3298 double vStari = vStar_dof[i];
3300 double u_beta_numerator = 0., u_beta_denominator = 0.;
3301 double v_beta_numerator = 0., v_beta_denominator = 0.;
3304 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3306 int j = colind_1D[offset];
3307 double uStarj = uStar_dof[j];
3308 double vStarj = vStar_dof[j];
3311 u_beta_numerator += (uStarj - uStari);
3312 u_beta_denominator += fabs(uStarj - uStari);
3314 v_beta_numerator += (vStarj - vStari);
3315 v_beta_denominator += fabs(vStarj - vStari);
3317 double u_beta = fabs(u_beta_numerator)/(u_beta_denominator+1E-10);
3318 double v_beta = fabs(v_beta_numerator)/(v_beta_denominator+1E-10);
3337 if (ARTIFICIAL_VISCOSITY==3)
3339 for(
int eN=0;eN<nElements_global;eN++)
3341 double uStar_He = element_uStar_He[eN];
3342 double vStar_He = element_vStar_He[eN];
3343 for(
int i=0;i<nDOF_test_element;i++)
3345 register int eN_i=eN*nDOF_test_element+i;
3346 register int gi = vel_l2g[eN_i];
3354 if (ARTIFICIAL_VISCOSITY==3)
3356 for (
int i=0; i<numDOFs_1D; i++)
3361 if (isBoundary_1D[i] == 1)
3385 for (
int i=0; i<numDOFs_1D; i++)
3388 double uStar_dii = 0;
3389 double vStar_dii = 0;
3390 double ui = u_dof[i];
3391 double vi = v_dof[i];
3393 double ith_u_dissipative_term = 0;
3394 double ith_v_dissipative_term = 0;
3399 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3401 int j = colind_1D[offset];
3404 double uj = u_dof[j];
3405 double vj = v_dof[j];
3410 if (ARTIFICIAL_VISCOSITY==4)
3412 double dEVij = fmax(laggedEntropyResidualPerNode[i],
3413 laggedEntropyResidualPerNode[j]);
3416 uStar_dMatrix[ij] = fmin(dLij,
cE*dEVij);
3417 vStar_dMatrix[i] = uStar_dMatrix[ij];
3426 uStar_dii -= uStar_dMatrix[ij];
3427 vStar_dii -= vStar_dMatrix[ij];
3429 ith_u_dissipative_term += uStar_dMatrix[ij]*(uj-ui);
3430 ith_v_dissipative_term += vStar_dMatrix[ij]*(vj-vi);
3439 uStar_dMatrix[ii] = uStar_dii;
3440 vStar_dMatrix[ii] = vStar_dii;
3441 globalResidual[offset_u+stride_u*i] += -ith_u_dissipative_term;
3442 globalResidual[offset_v+stride_v*i] += -ith_v_dissipative_term;
3453 std::memset(particle_netForces.data(),0,nParticles*3*
sizeof(
double));
3454 std::memset(particle_netMoments.data(),0,nParticles*3*
sizeof(
double));
3459 register double Fx = 0.0, Fy = 0.0, Fxp = 0.0, Fyp = 0.0, surfaceArea=0.0, Mz = 0.0;
3463 eN_nDOF_trial_element = eN*nDOF_trial_element;
3464 register double elementResidual_mesh[nDOF_test_element],
3465 elementResidual_p[nDOF_test_element],
3466 elementResidual_u[nDOF_test_element],
3467 elementResidual_v[nDOF_test_element],
3474 for (
int i=0;i<nDOF_test_element;i++)
3476 elementResidual_mesh[i]=0.0;
3477 elementResidual_p[i]=0.0;
3478 elementResidual_u[i]=0.0;
3479 elementResidual_v[i]=0.0;
3482 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3484 register int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
3486 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3487 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3495 jac_ext[nSpace*nSpace],
3497 jacInv_ext[nSpace*nSpace],
3498 boundaryJac[nSpace*(nSpace-1)],
3499 metricTensor[(nSpace-1)*(nSpace-1)],
3500 metricTensorDetSqrt,
3501 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
3502 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
3503 vel_grad_test_dS[nDOF_trial_element*nSpace],
3504 normal[2],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3505 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3506 force_x,force_y,force_z,force_p_x,force_p_y,force_p_z,force_v_x,force_v_y,force_v_z,r_x,r_y,r_z;
3508 ck.calculateMapping_elementBoundary(eN,
3514 mesh_trial_trace_ref.data(),
3515 mesh_grad_trial_trace_ref.data(),
3516 boundaryJac_ref.data(),
3522 metricTensorDetSqrt,
3526 ck.calculateMappingVelocity_elementBoundary(eN,
3530 mesh_velocity_dof.data(),
3532 mesh_trial_trace_ref.data(),
3533 xt_ext,yt_ext,zt_ext,
3538 dS = metricTensorDetSqrt*dS_ref[kb];
3540 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3543 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3545 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3546 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3548 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3549 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3551 for (
int j=0;j<nDOF_trial_element;j++)
3553 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3554 for (
int I=0;I<nSpace;I++)
3555 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3559 double distance[2], P_normal[2], P_tangent[2];
3563 if(use_ball_as_particle==1)
3571 P_normal[0],P_normal[1]);
3573 ball_velocity.data(),ball_angular_velocity.data(),
3575 x_ext-dist*P_normal[0],
3576 y_ext-dist*P_normal[1],
3582 dist = ebq_global_phi_solid[ebN_kb];
3583 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*3+0];
3584 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*3+1];
3585 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*3+0];
3586 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*3+1];
3590 ck.calculateGScale(G,normal,h_penalty);
3594 assert(h_penalty>0.0);
3595 if (h_penalty < std::abs(dist))
3596 h_penalty = std::abs(dist);
3597 distance[0] = -P_normal[0]*dist;
3598 distance[1] = -P_normal[1]*dist;
3599 P_tangent[0] = -P_normal[1];
3600 P_tangent[1] = P_normal[0];
3602 double C_adim =
C_sbm*visco/h_penalty;
3603 double beta_adim =
beta_sbm*visco/h_penalty;
3608 const double u_m_uD[2] = {u_ext - bc_u_ext,v_ext - bc_v_ext};
3609 const double zero_vec[2]={0.,0.};
3612 for (
int i=0;i<nDOF_test_element;i++)
3614 int eN_i = eN*nDOF_test_element+i;
3616 int GlobPos_u = offset_u+stride_u*vel_l2g[eN_i];
3617 int GlobPos_v = offset_v+stride_v*vel_l2g[eN_i];
3618 double phi_i = vel_test_dS[i];
3619 double Gxphi_i = vel_grad_test_dS[i*nSpace+0];
3620 double Gyphi_i = vel_grad_test_dS[i*nSpace+1];
3621 double *grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
3623 const double grad_phi_i_dot_t =
get_dot_product(P_tangent,grad_phi_i);
3626 globalResidual[GlobPos_u] += C_adim*phi_i*u_m_uD[0];
3627 globalResidual[GlobPos_v] += C_adim*phi_i*u_m_uD[1];
3628 Fx += C_adim*phi_i*u_m_uD[0];
3629 Fy += C_adim*phi_i*u_m_uD[1];
3633 globalResidual[GlobPos_u] -= visco * phi_i*res[0];
3634 globalResidual[GlobPos_v] -= visco * phi_i*res[1];
3635 Fx -= visco * phi_i*res[0];
3636 Fy -= visco * phi_i*res[1];
3649 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*u_m_uD[0];
3650 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*u_m_uD[1];
3651 Fx += C_adim*grad_phi_i_dot_d*u_m_uD[0];
3652 Fy += C_adim*grad_phi_i_dot_d*u_m_uD[1];
3655 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*grad_u_d[0];
3656 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*grad_u_d[1];
3657 Fx += C_adim*grad_phi_i_dot_d*grad_u_d[0];
3658 Fy += C_adim*grad_phi_i_dot_d*grad_u_d[1];
3661 globalResidual[GlobPos_u] += C_adim*phi_i*grad_u_d[0];
3662 globalResidual[GlobPos_v] += C_adim*phi_i*grad_u_d[1];
3663 Fx += C_adim*phi_i*grad_u_d[0];
3664 Fy += C_adim*phi_i*grad_u_d[1];
3678 globalResidual[GlobPos_u] += beta_adim*grad_u_t[0]*grad_phi_i_dot_t;
3679 globalResidual[GlobPos_v] += beta_adim*grad_u_t[1]*grad_phi_i_dot_t;
3680 Fx += beta_adim*grad_u_t[0]*grad_phi_i_dot_t;
3681 Fy += beta_adim*grad_u_t[1]*grad_phi_i_dot_t;
3690 for (
int i=0; i<nDOF_per_element_pressure;++i)
3692 p_ext += p_dof[p_l2g[eN*nDOF_per_element_pressure+i]]*p_trial_trace_ref[ebN_local_kb*nDOF_per_element_pressure+i];
3694 double nx = P_normal[0];
3695 double ny = P_normal[1];
3701 if(use_ball_as_particle==1)
3711 Mz += r_x*Fy-r_y*Fx;
3714 && ebN < nElementBoundaries_owned)
3735 gf_s.useExact=
false;
3736 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3738 register int ebN = exteriorElementBoundariesArray[ebNE],
3739 eN = elementBoundaryElementsArray[ebN*2+0],
3740 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
3741 eN_nDOF_trial_element = eN*nDOF_trial_element;
3742 register double elementResidual_mesh[nDOF_test_element],
3743 elementResidual_p[nDOF_test_element],
3744 elementResidual_u[nDOF_test_element],
3745 elementResidual_v[nDOF_test_element],
3748 const double* elementResidual_w(NULL);
3749 for (
int i=0;i<nDOF_test_element;i++)
3751 elementResidual_mesh[i]=0.0;
3752 elementResidual_p[i]=0.0;
3753 elementResidual_u[i]=0.0;
3754 elementResidual_v[i]=0.0;
3757 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3759 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3760 ebNE_kb_nSpace = ebNE_kb*nSpace,
3761 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3762 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3763 register double p_ext=0.0,
3772 dmom_u_acc_u_ext=0.0,
3774 dmom_v_acc_v_ext=0.0,
3776 dmom_w_acc_w_ext=0.0,
3777 mass_adv_ext[nSpace],
3778 dmass_adv_u_ext[nSpace],
3779 dmass_adv_v_ext[nSpace],
3780 dmass_adv_w_ext[nSpace],
3781 mom_u_adv_ext[nSpace],
3782 dmom_u_adv_u_ext[nSpace],
3783 dmom_u_adv_v_ext[nSpace],
3784 dmom_u_adv_w_ext[nSpace],
3785 mom_v_adv_ext[nSpace],
3786 dmom_v_adv_u_ext[nSpace],
3787 dmom_v_adv_v_ext[nSpace],
3788 dmom_v_adv_w_ext[nSpace],
3789 mom_w_adv_ext[nSpace],
3790 dmom_w_adv_u_ext[nSpace],
3791 dmom_w_adv_v_ext[nSpace],
3792 dmom_w_adv_w_ext[nSpace],
3793 mom_uu_diff_ten_ext[nSpace],
3794 mom_vv_diff_ten_ext[nSpace],
3795 mom_ww_diff_ten_ext[nSpace],
3796 mom_uv_diff_ten_ext[1],
3797 mom_uw_diff_ten_ext[1],
3798 mom_vu_diff_ten_ext[1],
3799 mom_vw_diff_ten_ext[1],
3800 mom_wu_diff_ten_ext[1],
3801 mom_wv_diff_ten_ext[1],
3802 mom_u_source_ext=0.0,
3803 mom_v_source_ext=0.0,
3804 mom_w_source_ext=0.0,
3806 dmom_u_ham_grad_p_ext[nSpace],
3807 dmom_u_ham_grad_u_ext[nSpace],
3809 dmom_v_ham_grad_p_ext[nSpace],
3810 dmom_v_ham_grad_v_ext[nSpace],
3812 dmom_w_ham_grad_p_ext[nSpace],
3813 dmom_w_ham_grad_w_ext[nSpace],
3814 dmom_u_adv_p_ext[nSpace],
3815 dmom_v_adv_p_ext[nSpace],
3816 dmom_w_adv_p_ext[nSpace],
3818 flux_mom_u_adv_ext=0.0,
3819 flux_mom_v_adv_ext=0.0,
3820 flux_mom_w_adv_ext=0.0,
3821 flux_mom_uu_diff_ext=0.0,
3822 flux_mom_uv_diff_ext=0.0,
3823 flux_mom_uw_diff_ext=0.0,
3824 flux_mom_vu_diff_ext=0.0,
3825 flux_mom_vv_diff_ext=0.0,
3826 flux_mom_vw_diff_ext=0.0,
3827 flux_mom_wu_diff_ext=0.0,
3828 flux_mom_wv_diff_ext=0.0,
3829 flux_mom_ww_diff_ext=0.0,
3834 bc_mom_u_acc_ext=0.0,
3835 bc_dmom_u_acc_u_ext=0.0,
3836 bc_mom_v_acc_ext=0.0,
3837 bc_dmom_v_acc_v_ext=0.0,
3838 bc_mom_w_acc_ext=0.0,
3839 bc_dmom_w_acc_w_ext=0.0,
3840 bc_mass_adv_ext[nSpace],
3841 bc_dmass_adv_u_ext[nSpace],
3842 bc_dmass_adv_v_ext[nSpace],
3843 bc_dmass_adv_w_ext[nSpace],
3844 bc_mom_u_adv_ext[nSpace],
3845 bc_dmom_u_adv_u_ext[nSpace],
3846 bc_dmom_u_adv_v_ext[nSpace],
3847 bc_dmom_u_adv_w_ext[nSpace],
3848 bc_mom_v_adv_ext[nSpace],
3849 bc_dmom_v_adv_u_ext[nSpace],
3850 bc_dmom_v_adv_v_ext[nSpace],
3851 bc_dmom_v_adv_w_ext[nSpace],
3852 bc_mom_w_adv_ext[nSpace],
3853 bc_dmom_w_adv_u_ext[nSpace],
3854 bc_dmom_w_adv_v_ext[nSpace],
3855 bc_dmom_w_adv_w_ext[nSpace],
3856 bc_mom_uu_diff_ten_ext[nSpace],
3857 bc_mom_vv_diff_ten_ext[nSpace],
3858 bc_mom_ww_diff_ten_ext[nSpace],
3859 bc_mom_uv_diff_ten_ext[1],
3860 bc_mom_uw_diff_ten_ext[1],
3861 bc_mom_vu_diff_ten_ext[1],
3862 bc_mom_vw_diff_ten_ext[1],
3863 bc_mom_wu_diff_ten_ext[1],
3864 bc_mom_wv_diff_ten_ext[1],
3865 bc_mom_u_source_ext=0.0,
3866 bc_mom_v_source_ext=0.0,
3867 bc_mom_w_source_ext=0.0,
3868 bc_mom_u_ham_ext=0.0,
3869 bc_dmom_u_ham_grad_p_ext[nSpace],
3870 bc_dmom_u_ham_grad_u_ext[nSpace],
3871 bc_mom_v_ham_ext=0.0,
3872 bc_dmom_v_ham_grad_p_ext[nSpace],
3873 bc_dmom_v_ham_grad_v_ext[nSpace],
3874 bc_mom_w_ham_ext=0.0,
3875 bc_dmom_w_ham_grad_p_ext[nSpace],
3876 bc_dmom_w_ham_grad_w_ext[nSpace],
3877 jac_ext[nSpace*nSpace],
3879 jacInv_ext[nSpace*nSpace],
3880 boundaryJac[nSpace*(nSpace-1)],
3881 metricTensor[(nSpace-1)*(nSpace-1)],
3882 metricTensorDetSqrt,
3883 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
3884 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
3885 vel_grad_test_dS[nDOF_trial_element*nSpace],
3886 normal[2],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3890 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3891 force_x,force_y,force_z,force_p_x,force_p_y,force_p_z,force_v_x,force_v_y,force_v_z,r_x,r_y,r_z;
3893 ck.calculateMapping_elementBoundary(eN,
3899 mesh_trial_trace_ref.data(),
3900 mesh_grad_trial_trace_ref.data(),
3901 boundaryJac_ref.data(),
3907 metricTensorDetSqrt,
3911 ck.calculateMappingVelocity_elementBoundary(eN,
3915 mesh_velocity_dof.data(),
3917 mesh_trial_trace_ref.data(),
3918 xt_ext,yt_ext,zt_ext,
3930 dS = metricTensorDetSqrt*dS_ref[kb];
3933 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3934 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
3936 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3937 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3938 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3943 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3947 p_ext = ebqe_p[ebNE_kb];
3948 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3949 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3952 for (
int I=0;I<nSpace;I++)
3953 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace + I];
3954 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3955 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3958 for (
int j=0;j<nDOF_trial_element;j++)
3961 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3962 for (
int I=0;I<nSpace;I++)
3963 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3965 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
3967 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
3968 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
3971 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
3975 double distance_to_omega_solid = 1e10;
3976 if (use_ball_as_particle == 1)
3978 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
3982 distance_to_omega_solid = ebq_global_phi_solid[ebN*nQuadraturePoints_elementBoundary+kb];
3984 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);
3993 elementDiameter[eN],
3994 smagorinskyConstant,
3995 turbulenceClosureModel,
3998 ebqe_vf_ext[ebNE_kb],
3999 ebqe_phi_ext[ebNE_kb],
4000 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
4001 distance_to_omega_solid,
4002 ebqe_kappa_phi_ext[ebNE_kb],
4014 ebqe_velocity_star[ebNE_kb_nSpace+0],
4015 ebqe_velocity_star[ebNE_kb_nSpace+1],
4016 ebqe_velocity_star[ebNE_kb_nSpace+1],
4040 mom_uu_diff_ten_ext,
4041 mom_vv_diff_ten_ext,
4042 mom_ww_diff_ten_ext,
4043 mom_uv_diff_ten_ext,
4044 mom_uw_diff_ten_ext,
4045 mom_vu_diff_ten_ext,
4046 mom_vw_diff_ten_ext,
4047 mom_wu_diff_ten_ext,
4048 mom_wv_diff_ten_ext,
4053 dmom_u_ham_grad_p_ext,
4054 dmom_u_ham_grad_u_ext,
4056 dmom_v_ham_grad_p_ext,
4057 dmom_v_ham_grad_v_ext,
4059 dmom_w_ham_grad_p_ext,
4060 dmom_w_ham_grad_w_ext,
4068 MATERIAL_PARAMETERS_AS_FUNCTION,
4069 ebqe_density_as_function[ebNE_kb],
4070 ebqe_dynamic_viscosity_as_function[ebNE_kb],
4073 use_ball_as_particle,
4076 ball_velocity.data(),
4077 ball_angular_velocity.data(),
4078 INT_BY_PARTS_PRESSURE);
4087 elementDiameter[eN],
4088 smagorinskyConstant,
4089 turbulenceClosureModel,
4092 bc_ebqe_vf_ext[ebNE_kb],
4093 bc_ebqe_phi_ext[ebNE_kb],
4094 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
4095 distance_to_omega_solid,
4096 ebqe_kappa_phi_ext[ebNE_kb],
4108 ebqe_velocity_star[ebNE_kb_nSpace+0],
4109 ebqe_velocity_star[ebNE_kb_nSpace+1],
4110 ebqe_velocity_star[ebNE_kb_nSpace+1],
4111 bc_eddy_viscosity_ext,
4113 bc_dmom_u_acc_u_ext,
4115 bc_dmom_v_acc_v_ext,
4117 bc_dmom_w_acc_w_ext,
4123 bc_dmom_u_adv_u_ext,
4124 bc_dmom_u_adv_v_ext,
4125 bc_dmom_u_adv_w_ext,
4127 bc_dmom_v_adv_u_ext,
4128 bc_dmom_v_adv_v_ext,
4129 bc_dmom_v_adv_w_ext,
4131 bc_dmom_w_adv_u_ext,
4132 bc_dmom_w_adv_v_ext,
4133 bc_dmom_w_adv_w_ext,
4134 bc_mom_uu_diff_ten_ext,
4135 bc_mom_vv_diff_ten_ext,
4136 bc_mom_ww_diff_ten_ext,
4137 bc_mom_uv_diff_ten_ext,
4138 bc_mom_uw_diff_ten_ext,
4139 bc_mom_vu_diff_ten_ext,
4140 bc_mom_vw_diff_ten_ext,
4141 bc_mom_wu_diff_ten_ext,
4142 bc_mom_wv_diff_ten_ext,
4143 bc_mom_u_source_ext,
4144 bc_mom_v_source_ext,
4145 bc_mom_w_source_ext,
4147 bc_dmom_u_ham_grad_p_ext,
4148 bc_dmom_u_ham_grad_u_ext,
4150 bc_dmom_v_ham_grad_p_ext,
4151 bc_dmom_v_ham_grad_v_ext,
4153 bc_dmom_w_ham_grad_p_ext,
4154 bc_dmom_w_ham_grad_w_ext,
4162 MATERIAL_PARAMETERS_AS_FUNCTION,
4163 ebqe_density_as_function[ebNE_kb],
4164 ebqe_dynamic_viscosity_as_function[ebNE_kb],
4167 use_ball_as_particle,
4170 ball_velocity.data(),
4171 ball_angular_velocity.data(),
4172 INT_BY_PARTS_PRESSURE);
4175 if (turbulenceClosureModel >= 3)
4177 const double turb_var_grad_0_dummy[2] = {0.,0.};
4178 const double c_mu = 0.09;
4187 ebqe_vf_ext[ebNE_kb],
4188 ebqe_phi_ext[ebNE_kb],
4191 ebqe_turb_var_0[ebNE_kb],
4192 ebqe_turb_var_1[ebNE_kb],
4193 turb_var_grad_0_dummy,
4195 mom_uu_diff_ten_ext,
4196 mom_vv_diff_ten_ext,
4197 mom_ww_diff_ten_ext,
4198 mom_uv_diff_ten_ext,
4199 mom_uw_diff_ten_ext,
4200 mom_vu_diff_ten_ext,
4201 mom_vw_diff_ten_ext,
4202 mom_wu_diff_ten_ext,
4203 mom_wv_diff_ten_ext,
4216 bc_ebqe_vf_ext[ebNE_kb],
4217 bc_ebqe_phi_ext[ebNE_kb],
4220 ebqe_turb_var_0[ebNE_kb],
4221 ebqe_turb_var_1[ebNE_kb],
4222 turb_var_grad_0_dummy,
4223 bc_eddy_viscosity_ext,
4224 bc_mom_uu_diff_ten_ext,
4225 bc_mom_vv_diff_ten_ext,
4226 bc_mom_ww_diff_ten_ext,
4227 bc_mom_uv_diff_ten_ext,
4228 bc_mom_uw_diff_ten_ext,
4229 bc_mom_vu_diff_ten_ext,
4230 bc_mom_vw_diff_ten_ext,
4231 bc_mom_wu_diff_ten_ext,
4232 bc_mom_wv_diff_ten_ext,
4233 bc_mom_u_source_ext,
4234 bc_mom_v_source_ext,
4235 bc_mom_w_source_ext);
4242 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
4243 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
4245 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
4246 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
4249 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
4250 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
4252 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
4253 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
4265 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
4266 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
4269 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
4270 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
4279 ck.calculateGScale(G,normal,h_penalty);
4280 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
4282 isDOFBoundary_u[ebNE_kb],
4283 isDOFBoundary_v[ebNE_kb],
4284 isDOFBoundary_w[ebNE_kb],
4285 isAdvectiveFluxBoundary_p[ebNE_kb],
4286 isAdvectiveFluxBoundary_u[ebNE_kb],
4287 isAdvectiveFluxBoundary_v[ebNE_kb],
4288 isAdvectiveFluxBoundary_w[ebNE_kb],
4289 dmom_u_ham_grad_p_ext[0],
4290 bc_dmom_u_ham_grad_p_ext[0],
4301 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),
4302 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
4303 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
4304 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
4332 &ebqe_velocity_star[ebNE_kb_nSpace],
4333 &ebqe_velocity[ebNE_kb_nSpace]);
4335 for (
int I=0;I<nSpace;I++)
4337 ebqe_grad_u[ebNE_kb_nSpace+I] = grad_u_ext[I];
4338 ebqe_grad_v[ebNE_kb_nSpace+I] = grad_v_ext[I];
4342 ebqe_phi_ext[ebNE_kb],
4343 sdInfo_u_u_rowptr.data(),
4344 sdInfo_u_u_colind.data(),
4345 isDOFBoundary_u[ebNE_kb],
4346 isDiffusiveFluxBoundary_u[ebNE_kb],
4348 bc_mom_uu_diff_ten_ext,
4350 ebqe_bc_flux_u_diff_ext[ebNE_kb],
4351 mom_uu_diff_ten_ext,
4355 flux_mom_uu_diff_ext);
4357 ebqe_phi_ext[ebNE_kb],
4358 sdInfo_u_v_rowptr.data(),
4359 sdInfo_u_v_colind.data(),
4360 isDOFBoundary_v[ebNE_kb],
4361 isDiffusiveFluxBoundary_v[ebNE_kb],
4363 bc_mom_uv_diff_ten_ext,
4366 mom_uv_diff_ten_ext,
4370 flux_mom_uv_diff_ext);
4387 ebqe_phi_ext[ebNE_kb],
4388 sdInfo_v_u_rowptr.data(),
4389 sdInfo_v_u_colind.data(),
4390 isDOFBoundary_u[ebNE_kb],
4391 isDiffusiveFluxBoundary_u[ebNE_kb],
4393 bc_mom_vu_diff_ten_ext,
4396 mom_vu_diff_ten_ext,
4400 flux_mom_vu_diff_ext);
4402 ebqe_phi_ext[ebNE_kb],
4403 sdInfo_v_v_rowptr.data(),
4404 sdInfo_v_v_colind.data(),
4405 isDOFBoundary_v[ebNE_kb],
4406 isDiffusiveFluxBoundary_v[ebNE_kb],
4408 bc_mom_vv_diff_ten_ext,
4410 ebqe_bc_flux_v_diff_ext[ebNE_kb],
4411 mom_vv_diff_ten_ext,
4415 flux_mom_vv_diff_ext);
4476 flux[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
4484 if (ebN < nElementBoundaries_owned)
4486 force_v_x = (flux_mom_u_adv_ext + flux_mom_uu_diff_ext + flux_mom_uv_diff_ext + flux_mom_uw_diff_ext)/dmom_u_ham_grad_p_ext[0];
4487 force_v_y = (flux_mom_v_adv_ext + flux_mom_vu_diff_ext + flux_mom_vv_diff_ext + flux_mom_vw_diff_ext)/dmom_u_ham_grad_p_ext[0];
4490 force_p_x = p_ext*normal[0];
4491 force_p_y = p_ext*normal[1];
4494 force_x = force_p_x + force_v_x;
4495 force_y = force_p_y + force_v_y;
4498 r_x = x_ext - barycenters[3*boundaryFlags[ebN]+0];
4499 r_y = y_ext - barycenters[3*boundaryFlags[ebN]+1];
4502 wettedAreas[boundaryFlags[ebN]] += dS*(1.0-ebqe_vf_ext[ebNE_kb]);
4504 netForces_p[3*boundaryFlags[ebN]+0] += force_p_x*dS;
4505 netForces_p[3*boundaryFlags[ebN]+1] += force_p_y*dS;
4508 netForces_v[3*boundaryFlags[ebN]+0] += force_v_x*dS;
4509 netForces_v[3*boundaryFlags[ebN]+1] += force_v_y*dS;
4514 netMoments[3*boundaryFlags[ebN]+2] += (r_x*force_y - r_y*force_x)*dS;
4519 for (
int i=0;i<nDOF_test_element;i++)
4525 elementResidual_u[i] +=
4526 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[0] : 0.) +
4527 ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
4528 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
4529 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
4530 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4531 isDiffusiveFluxBoundary_u[ebNE_kb],
4536 sdInfo_u_u_rowptr.data(),
4537 sdInfo_u_u_colind.data(),
4538 mom_uu_diff_ten_ext,
4539 &vel_grad_test_dS[i*nSpace])+
4540 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4541 isDiffusiveFluxBoundary_u[ebNE_kb],
4546 sdInfo_u_v_rowptr.data(),
4547 sdInfo_u_v_colind.data(),
4548 mom_uv_diff_ten_ext,
4549 &vel_grad_test_dS[i*nSpace]);
4560 elementResidual_v[i] +=
4561 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[1] : 0.) +
4562 ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
4563 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
4564 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
4565 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4566 isDiffusiveFluxBoundary_v[ebNE_kb],
4571 sdInfo_v_u_rowptr.data(),
4572 sdInfo_v_u_colind.data(),
4573 mom_vu_diff_ten_ext,
4574 &vel_grad_test_dS[i*nSpace])+
4575 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4576 isDiffusiveFluxBoundary_v[ebNE_kb],
4581 sdInfo_v_v_rowptr.data(),
4582 sdInfo_v_v_colind.data(),
4583 mom_vv_diff_ten_ext,
4584 &vel_grad_test_dS[i*nSpace]);
4637 for (
int i=0;i<nDOF_test_element;i++)
4639 int eN_i = eN*nDOF_test_element+i;
4644 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=elementResidual_u[i];
4645 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=elementResidual_v[i];
4649 gf.useExact=useExact;
4650 gf_s.useExact=useExact;
4657 particle_surfaceArea[0] = cut_cell_boundary_length;
4658 particle_netForces[(0+nParticles)*3 +0] = p_force_x;
4659 particle_netForces[(0+nParticles)*3 +1] = p_force_y;
4660 std::cout<<
"===end mesh==="<<std::endl<<std::flush;
4667 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4668 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4669 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4670 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
4671 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
4672 double PSTAB = args.
scalar<
double>(
"PSTAB");
4673 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4674 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
4675 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4676 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4677 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4678 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4679 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4680 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
4681 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
4682 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
4683 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
4684 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4685 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4686 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
4687 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4688 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4689 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4690 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4691 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4692 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4693 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4694 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4695 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4696 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4697 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4698 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4699 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4700 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4701 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4702 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
4703 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4704 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4705 double hFactor = args.
scalar<
double>(
"hFactor");
4706 int nElements_global = args.
scalar<
int>(
"nElements_global");
4707 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
4708 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
4709 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
4710 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
4711 double useRBLES = args.
scalar<
double>(
"useRBLES");
4712 double useMetrics = args.
scalar<
double>(
"useMetrics");
4713 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4714 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4715 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4716 double sigma = args.
scalar<
double>(
"sigma");
4721 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4722 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4723 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4724 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4725 double C_dg = args.
scalar<
double>(
"C_dg");
4726 double C_b = args.
scalar<
double>(
"C_b");
4727 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
4728 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
4729 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
4730 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
4731 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
4732 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
4733 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
4734 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
4735 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
4736 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
4737 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
4738 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
4739 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
4740 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
4741 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
4742 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
4743 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
4744 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4745 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4746 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4747 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4748 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4749 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4750 xt::pyarray<double>& g = args.
array<
double>(
"g");
4751 const double useVF = args.
scalar<
double>(
"useVF");
4752 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4753 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4754 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
4755 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4756 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4757 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4758 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4759 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4760 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
4761 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
4762 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4763 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
4764 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4765 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4766 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4767 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4768 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4769 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4770 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4771 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4772 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4773 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4774 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4775 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4776 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4777 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4778 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4779 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4780 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4781 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4782 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4783 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4784 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4785 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4786 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
4787 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
4788 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
4789 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
4790 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
4791 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
4792 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
4793 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
4794 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
4795 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
4796 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
4797 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
4798 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
4799 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
4800 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
4801 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
4802 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
4803 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
4804 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
4805 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
4806 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
4807 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
4808 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
4809 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
4810 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
4811 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
4812 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
4813 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
4814 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
4815 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
4816 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
4817 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
4818 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
4819 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4820 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4821 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
4822 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4823 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4824 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4825 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
4826 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4827 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
4828 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4829 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4830 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
4831 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
4832 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
4833 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4834 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4835 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4836 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4837 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4838 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4839 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4840 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4841 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4842 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4843 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4844 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4845 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4846 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4847 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4848 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4849 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4850 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4851 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4852 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4853 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4854 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4855 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4856 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
4857 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
4858 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
4859 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
4860 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
4861 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
4862 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
4863 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
4864 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
4865 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
4866 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
4867 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
4868 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
4869 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
4870 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
4871 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
4872 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
4873 int nParticles = args.
scalar<
int>(
"nParticles");
4874 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
4875 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
4876 double particle_beta = args.
scalar<
double>(
"particle_beta");
4877 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
4878 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
4879 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
4880 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
4881 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
4882 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
4883 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
4884 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
4885 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
4886 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
4887 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
4888 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
4889 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
4890 double dt = args.
scalar<
double>(
"dt");
4891 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
4892 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
4893 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
4894 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
4895 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
4896 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
4897 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
4898 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
4899 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
4900 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
4901 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
4902 int offset_u = args.
scalar<
int>(
"offset_u");
4903 int offset_v = args.
scalar<
int>(
"offset_v");
4904 int offset_w = args.
scalar<
int>(
"offset_w");
4905 int stride_u = args.
scalar<
int>(
"stride_u");
4906 int stride_v = args.
scalar<
int>(
"stride_v");
4907 int stride_w = args.
scalar<
int>(
"stride_w");
4908 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
4909 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
4910 xt::pyarray<int>& rowptr = args.
array<
int>(
"rowptr");
4911 xt::pyarray<int>& colind = args.
array<
int>(
"colind");
4912 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
4916 std::valarray<double> particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3);
4917 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
4918 for(
int eN=0;eN<nElements_global;eN++)
4920 register double eps_rho,eps_mu;
4921 double element_active=1.0;
4923 register double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
4924 elementJacobian_p_u[nDOF_test_element][nDOF_trial_element],
4925 elementJacobian_p_v[nDOF_test_element][nDOF_trial_element],
4926 elementJacobian_p_w[nDOF_test_element][nDOF_trial_element],
4927 elementJacobian_u_p[nDOF_test_element][nDOF_trial_element],
4928 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
4929 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
4930 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
4931 elementJacobian_v_p[nDOF_test_element][nDOF_trial_element],
4932 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
4933 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
4934 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
4935 elementJacobian_w_p[nDOF_test_element][nDOF_trial_element],
4936 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
4937 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
4938 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
4939 for (
int i=0;i<nDOF_test_element;i++)
4940 for (
int j=0;j<nDOF_trial_element;j++)
4942 elementJacobian_p_p[i][j]=0.0;
4943 elementJacobian_p_u[i][j]=0.0;
4944 elementJacobian_p_v[i][j]=0.0;
4945 elementJacobian_p_w[i][j]=0.0;
4946 elementJacobian_u_p[i][j]=0.0;
4947 elementJacobian_u_u[i][j]=0.0;
4948 elementJacobian_u_v[i][j]=0.0;
4949 elementJacobian_u_w[i][j]=0.0;
4950 elementJacobian_v_p[i][j]=0.0;
4951 elementJacobian_v_u[i][j]=0.0;
4952 elementJacobian_v_v[i][j]=0.0;
4953 elementJacobian_v_w[i][j]=0.0;
4954 elementJacobian_w_p[i][j]=0.0;
4955 elementJacobian_w_u[i][j]=0.0;
4956 elementJacobian_w_v[i][j]=0.0;
4957 elementJacobian_w_w[i][j]=0.0;
4968 double _distance[nDOF_mesh_trial_element]={0.0};
4970 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4972 if(use_ball_as_particle==1)
4975 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
4976 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
4977 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
4982 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
4984 if ( _distance[I] >= 0)
4987 if (pos_counter == 2)
4992 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4994 if (_distance[I] < 0)
4997 assert(opp_node >=0);
4998 assert(opp_node <nDOF_mesh_trial_element);
5000 else if (pos_counter == 3)
5009 if(use_ball_as_particle==1)
5011 for (
int I=0;I<nDOF_mesh_trial_element;I++)
5013 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
5014 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
5015 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
5016 phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]]);
5018 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
5019 for (
int j=0;j<nDOF_mesh_trial_element;j++)
5021 register int eN_j = eN*nDOF_mesh_trial_element+j;
5022 element_phi[j] = phi_dof[p_l2g[eN_j]];
5023 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
5025 double element_nodes[nDOF_mesh_trial_element*3];
5026 for (
int i=0;i<nDOF_mesh_trial_element;i++)
5028 register int eN_i=eN*nDOF_mesh_trial_element+i;
5029 for(
int I=0;I<3;I++)
5030 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
5032 gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),
false);
5033 gf.calculate(element_phi, element_nodes, x_ref.data(),
false);
5034 for (
int k=0;k<nQuadraturePoints_element;k++)
5038 int eN_k = eN*nQuadraturePoints_element+k,
5039 eN_k_nSpace = eN_k*nSpace,
5041 eN_nDOF_trial_element = eN*nDOF_trial_element;
5044 register double p=0.0,
u=0.0,
v=0.0,
w=0.0,
5045 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
5054 dmass_adv_u[nSpace],
5055 dmass_adv_v[nSpace],
5056 dmass_adv_w[nSpace],
5058 dmom_u_adv_u[nSpace],
5059 dmom_u_adv_v[nSpace],
5060 dmom_u_adv_w[nSpace],
5062 dmom_v_adv_u[nSpace],
5063 dmom_v_adv_v[nSpace],
5064 dmom_v_adv_w[nSpace],
5066 dmom_w_adv_u[nSpace],
5067 dmom_w_adv_v[nSpace],
5068 dmom_w_adv_w[nSpace],
5069 mom_uu_diff_ten[nSpace],
5070 mom_vv_diff_ten[nSpace],
5071 mom_ww_diff_ten[nSpace],
5082 dmom_u_ham_grad_p[nSpace],
5083 dmom_u_ham_grad_u[nSpace],
5085 dmom_v_ham_grad_p[nSpace],
5086 dmom_v_ham_grad_v[nSpace],
5088 dmom_w_ham_grad_p[nSpace],
5089 dmom_w_ham_grad_w[nSpace],
5100 dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
5101 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
5102 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
5103 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
5104 Lstar_u_p[nDOF_test_element],
5105 Lstar_v_p[nDOF_test_element],
5106 Lstar_w_p[nDOF_test_element],
5107 Lstar_u_u[nDOF_test_element],
5108 Lstar_v_v[nDOF_test_element],
5109 Lstar_w_w[nDOF_test_element],
5110 Lstar_p_u[nDOF_test_element],
5111 Lstar_p_v[nDOF_test_element],
5112 Lstar_p_w[nDOF_test_element],
5117 dsubgridError_p_u[nDOF_trial_element],
5118 dsubgridError_p_v[nDOF_trial_element],
5119 dsubgridError_p_w[nDOF_trial_element],
5120 dsubgridError_u_p[nDOF_trial_element],
5121 dsubgridError_u_u[nDOF_trial_element],
5122 dsubgridError_v_p[nDOF_trial_element],
5123 dsubgridError_v_v[nDOF_trial_element],
5124 dsubgridError_w_p[nDOF_trial_element],
5125 dsubgridError_w_w[nDOF_trial_element],
5126 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
5127 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
5130 jacInv[nSpace*nSpace],
5131 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
5132 vel_hess_trial[nDOF_trial_element*
nSpace2],
5134 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_test_element],
5135 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
5140 dmom_u_source[nSpace],
5141 dmom_v_source[nSpace],
5142 dmom_w_source[nSpace],
5145 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
5147 ck.calculateMapping_element(eN,
5151 mesh_trial_ref.data(),
5152 mesh_grad_trial_ref.data(),
5157 ck.calculateH_element(eN,
5159 nodeDiametersArray.data(),
5161 mesh_trial_ref.data(),
5163 ck.calculateMappingVelocity_element(eN,
5165 mesh_velocity_dof.data(),
5167 mesh_trial_ref.data(),
5172 dV = fabs(jacDet)*dV_ref[k];
5173 ck.calculateG(jacInv,G,G_dd_G,tr_G);
5176 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
5177 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
5178 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
5182 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
5183 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
5187 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
5188 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
5192 for (
int I=0;I<nSpace;I++)
5193 grad_p[I] = q_grad_p[eN_k_nSpace+I];
5194 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
5195 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
5196 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
5197 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
5200 for (
int j=0;j<nDOF_trial_element;j++)
5203 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
5204 for (
int I=0;I<nSpace;I++)
5207 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
5211 double div_mesh_velocity=0.0;
5212 int NDOF_MESH_TRIAL_ELEMENT=3;
5213 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
5215 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
5216 div_mesh_velocity +=
5217 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*2+0] +
5218 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*2+1];
5220 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
5223 porosity = 1.0 - q_vos[eN_k];
5228 double distance_to_omega_solid = phi_solid[eN_k];
5229 double eddy_viscosity(0.),rhoSave,nuSave;
5238 elementDiameter[eN],
5239 smagorinskyConstant,
5240 turbulenceClosureModel,
5245 &normal_phi[eN_k_nSpace],
5246 distance_to_omega_solid,
5259 q_velocity_sge[eN_k_nSpace+0],
5260 q_velocity_sge[eN_k_nSpace+1],
5261 q_velocity_sge[eN_k_nSpace+1],
5313 MATERIAL_PARAMETERS_AS_FUNCTION,
5314 density_as_function[eN_k],
5315 dynamic_viscosity_as_function[eN_k],
5318 use_ball_as_particle,
5321 ball_velocity.data(),
5322 ball_angular_velocity.data(),
5323 INT_BY_PARTS_PRESSURE);
5325 mass_source = q_mass_source[eN_k];
5326 for (
int I=0;I<nSpace;I++)
5328 dmom_u_source[I] = 0.0;
5329 dmom_v_source[I] = 0.0;
5330 dmom_w_source[I] = 0.0;
5351 q_velocity_sge[eN_k_nSpace+0],
5352 q_velocity_sge[eN_k_nSpace+1],
5353 q_velocity_sge[eN_k_nSpace+1],
5354 eps_solid[elementFlags[eN]],
5356 q_velocity_solid[eN_k_nSpace+0],
5357 q_velocity_solid[eN_k_nSpace+1],
5358 q_velocity_solid[eN_k_nSpace+1],
5359 q_velocityStar_solid[eN_k_nSpace+0],
5360 q_velocityStar_solid[eN_k_nSpace+1],
5361 q_velocityStar_solid[eN_k_nSpace+1],
5368 q_grad_vos[eN_k_nSpace+0],
5369 q_grad_vos[eN_k_nSpace+1],
5370 q_grad_vos[eN_k_nSpace+1]);
5372 double C_particles=0.0;
5373 if(nParticles > 0 && USE_SBM==0)
5378 nQuadraturePoints_global,
5379 &particle_signed_distances[eN_k],
5380 &particle_signed_distance_normals[eN_k_3d],
5381 &particle_velocities[eN_k_3d],
5382 particle_centroids.data(),
5383 use_ball_as_particle,
5386 ball_velocity.data(),
5387 ball_angular_velocity.data(),
5389 particle_penalty_constant/h_phi,
5390 particle_alpha/h_phi,
5391 particle_beta/h_phi,
5408 q_velocity_sge[eN_k_nSpace+0],
5409 q_velocity_sge[eN_k_nSpace+1],
5410 q_velocity_sge[eN_k_nSpace+1],
5433 &particle_netForces[0],
5434 &particle_netMoments[0],
5435 &particle_surfaceArea[0]);
5437 if (turbulenceClosureModel >= 3)
5439 const double c_mu = 0.09;
5454 &q_turb_var_grad_0[eN_k_nSpace],
5474 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
5475 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
5477 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
5478 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
5481 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
5482 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
5484 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
5485 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
5498 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5504 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5518 mom_u_acc_t *= dmom_u_acc_u;
5519 mom_v_acc_t *= dmom_v_acc_v;
5522 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
5523 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
5529 ck.Mass_strong(-q_dvos_dt[eN_k]) +
5530 ck.Advection_strong(dmass_adv_u,grad_u) +
5531 ck.Advection_strong(dmass_adv_v,grad_v) +
5533 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
5535 ck.Reaction_strong(mass_source);
5539 ck.Mass_strong(mom_u_acc_t) +
5540 ck.Advection_strong(dmom_adv_sge,grad_u) +
5541 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
5542 ck.Reaction_strong(mom_u_source) -
5543 ck.Reaction_strong(
u*div_mesh_velocity);
5546 ck.Mass_strong(mom_v_acc_t) +
5547 ck.Advection_strong(dmom_adv_sge,grad_v) +
5548 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
5549 ck.Reaction_strong(mom_v_source) -
5550 ck.Reaction_strong(
v*div_mesh_velocity);
5560 for (
int j=0;j<nDOF_trial_element;j++)
5562 register int j_nSpace = j*nSpace;
5563 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
5564 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
5567 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
5568 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5569 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5570 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5572 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
5573 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5574 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5575 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5583 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j]);
5584 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j]);
5590 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
5592 elementDiameter[eN],
5597 dmom_u_ham_grad_p[0],
5607 dmom_u_ham_grad_p[0],
5613 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
5614 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
5647 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
5648 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
5652 for (
int i=0;i<nDOF_test_element;i++)
5654 register int i_nSpace = i*nSpace;
5655 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
5656 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
5658 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5659 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5661 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5662 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5665 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
5666 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
5671 dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5672 dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5675 dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5676 dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5684 double unit_normal[nSpace];
5685 double norm_grad_phi = 0.;
5686 for (
int I=0;I<nSpace;I++)
5687 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
5688 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
5689 for (
int I=0;I<nSpace;I++)
5690 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
5691 double delta =
gf.D(eps_mu,
phi[eN_k]);
5692 register double vel_tgrad_test_i[nSpace], vel_tgrad_test_j[nSpace];
5696 for(
int i=0;i<nDOF_test_element;i++)
5698 register int i_nSpace = i*nSpace;
5700 &vel_grad_trial[i_nSpace],
5702 for(
int j=0;j<nDOF_trial_element;j++)
5704 register int j_nSpace = j*nSpace;
5706 &vel_grad_trial[j_nSpace],
5722 elementJacobian_u_u[i][j] +=
5723 ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5724 ck.HamiltonianJacobian_weak(dmom_u_ham_grad_u,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5725 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5726 ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5728 ck.ReactionJacobian_weak(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5731 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
5732 ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5734 ck.NumericalDiffusion(dt*delta*sigma*dV,
5738 elementJacobian_u_v[i][j] +=
5739 ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5740 ck.SimpleDiffusionJacobian_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5742 ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5754 elementJacobian_v_u[i][j] +=
5755 ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5756 ck.SimpleDiffusionJacobian_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5758 ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5761 elementJacobian_v_v[i][j] +=
5762 ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5763 ck.HamiltonianJacobian_weak(dmom_v_ham_grad_v,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5764 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5765 ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5767 ck.ReactionJacobian_weak(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5770 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
5771 ck.NumericalDiffusionJacobian(q_numDiff_v_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5773 ck.NumericalDiffusion(dt*delta*sigma*dV,
5814 for (
int i=0;i<nDOF_test_element;i++)
5816 register int eN_i = eN*nDOF_test_element+i;
5817 for (
int j=0;j<nDOF_trial_element;j++)
5819 register int eN_i_j = eN_i*nDOF_trial_element+j;
5826 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += element_active*elementJacobian_u_u[i][j];
5827 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_u_v[eN_i_j]] += element_active*elementJacobian_u_v[i][j];
5831 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_v_u[eN_i_j]] += element_active*elementJacobian_v_u[i][j];
5832 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_v_v[eN_i_j]] += element_active*elementJacobian_v_v[i][j];
5844 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
5847 for (
int i=0; i<numDOFs_1D; i++)
5850 int u_gi = offset_u+stride_u*i;
5851 int v_gi = offset_v+stride_v*i;
5854 int u_ith_row_ptr = rowptr[u_gi];
5855 int v_ith_row_ptr = rowptr[v_gi];
5858 int numDOFs_ith_row = rowptr_1D[i+1]-rowptr_1D[i];
5859 for (
int counter = 0; counter < numDOFs_ith_row; counter++)
5862 int uu_ij = u_ith_row_ptr + (offset_u + counter*stride_u);
5863 int vv_ij = v_ith_row_ptr + (offset_v + counter*stride_v);
5866 double uStar_dij = uStar_dMatrix[ij];
5867 double vStar_dij = vStar_dMatrix[ij];
5870 globalJacobian[uu_ij] -= uStar_dij;
5871 globalJacobian[vv_ij] -= vStar_dij;
5888 eN_nDOF_trial_element = eN*nDOF_trial_element;
5889 register double eps_rho,eps_mu;
5893 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5895 register int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
5896 ebN_kb_nSpace = ebN_kb*nSpace,
5897 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5898 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5900 register double u_ext=0.0,
5906 jac_ext[nSpace*nSpace],
5908 jacInv_ext[nSpace*nSpace],
5909 boundaryJac[nSpace*(nSpace-1)],
5910 metricTensor[(nSpace-1)*(nSpace-1)],
5911 metricTensorDetSqrt,
5912 vel_grad_trial_trace[nDOF_trial_element*nSpace],
5914 vel_test_dS[nDOF_test_element],
5916 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
5917 vel_grad_test_dS[nDOF_trial_element*nSpace],
5918 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
5919 ck.calculateMapping_elementBoundary(eN,
5925 mesh_trial_trace_ref.data(),
5926 mesh_grad_trial_trace_ref.data(),
5927 boundaryJac_ref.data(),
5933 metricTensorDetSqrt,
5937 ck.calculateMappingVelocity_elementBoundary(eN,
5941 mesh_velocity_dof.data(),
5943 mesh_trial_trace_ref.data(),
5944 xt_ext,yt_ext,zt_ext,
5949 dS = metricTensorDetSqrt*dS_ref[kb];
5950 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
5953 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
5955 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
5956 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
5958 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
5959 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
5961 for (
int j=0;j<nDOF_trial_element;j++)
5963 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
5964 for (
int I=0;I<nSpace;I++)
5965 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
5972 ck.calculateGScale(G,normal,h_penalty);
5978 double distance[2], P_normal[2], P_tangent[2];
5980 if(use_ball_as_particle==1)
5988 P_normal[0],P_normal[1]);
5990 ball_velocity.data(), ball_angular_velocity.data(),
5992 x_ext-dist*P_normal[0],
5993 y_ext-dist*P_normal[1],
5999 dist = ebq_global_phi_solid[ebN_kb];
6000 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*3+0];
6001 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*3+1];
6002 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*3+0];
6003 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*3+1];
6005 distance[0] = -P_normal[0]*dist;
6006 distance[1] = -P_normal[1]*dist;
6007 P_tangent[0]= -P_normal[1];
6008 P_tangent[1]= P_normal[0];
6009 assert(h_penalty>0.0);
6010 if (h_penalty < std::abs(dist))
6011 h_penalty = std::abs(dist);
6014 double C_adim =
C_sbm*visco/h_penalty;
6015 double beta_adim =
beta_sbm*visco/h_penalty;
6017 for (
int i=0;i<nDOF_test_element;i++)
6019 register int eN_i = eN*nDOF_test_element+i;
6020 double phi_i = vel_test_dS[i];
6021 double* grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
6023 const double grad_phi_i_dot_t =
get_dot_product(P_tangent,grad_phi_i);
6026 const double zero_vec[2]={0.,0.};
6027 for (
int j=0;j<nDOF_trial_element;j++)
6032 + i*nDOF_trial_element
6035 double phi_j = vel_test_dS[j]/dS;
6036 const double grad_phi_j[2]={vel_grad_test_dS[j*nSpace+0]/dS,
6037 vel_grad_test_dS[j*nSpace+1]/dS};
6038 const double grad_phi_j_dot_d =
get_dot_product(distance, grad_phi_j);
6039 const double grad_phi_j_dot_t =
get_dot_product(P_tangent,grad_phi_j);
6043 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6045 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6050 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
6051 visco * phi_i * res[0];
6052 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
6053 visco * phi_i * res[1];
6056 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
6057 visco * phi_i * res[0];
6058 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
6059 visco * phi_i * res[1];
6063 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
6064 visco * phi_j * res[0];
6065 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
6066 visco * phi_j * res[1];
6068 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
6069 visco * phi_j * res[0];
6070 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
6071 visco * phi_j * res[1];
6074 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6075 C_adim*grad_phi_i_dot_d*phi_j;
6076 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6077 C_adim*grad_phi_i_dot_d*phi_j;
6080 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6081 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
6082 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6083 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
6086 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6087 C_adim*grad_phi_j_dot_d*phi_i;
6088 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6089 C_adim*grad_phi_j_dot_d*phi_i;
6093 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
6094 visco * grad_phi_j_dot_d * res[0];
6095 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
6096 visco * grad_phi_j_dot_d * res[1];
6099 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
6100 visco * grad_phi_j_dot_d * res[0] ;
6101 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
6102 visco * grad_phi_j_dot_d * res[1];
6107 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6108 beta_adim*grad_phi_j_dot_t*grad_phi_i_dot_t;
6109 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6110 beta_adim*grad_phi_j_dot_t*grad_phi_i_dot_t;
6122 gf_s.useExact=
false;
6123 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
6125 register int ebN = exteriorElementBoundariesArray[ebNE],
6126 eN = elementBoundaryElementsArray[ebN*2+0],
6127 eN_nDOF_trial_element = eN*nDOF_trial_element,
6128 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0];
6129 register double eps_rho,eps_mu;
6130 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
6132 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
6133 ebNE_kb_nSpace = ebNE_kb*nSpace,
6134 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
6135 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
6137 register double p_ext=0.0,
6146 dmom_u_acc_u_ext=0.0,
6148 dmom_v_acc_v_ext=0.0,
6150 dmom_w_acc_w_ext=0.0,
6151 mass_adv_ext[nSpace],
6152 dmass_adv_u_ext[nSpace],
6153 dmass_adv_v_ext[nSpace],
6154 dmass_adv_w_ext[nSpace],
6155 mom_u_adv_ext[nSpace],
6156 dmom_u_adv_u_ext[nSpace],
6157 dmom_u_adv_v_ext[nSpace],
6158 dmom_u_adv_w_ext[nSpace],
6159 mom_v_adv_ext[nSpace],
6160 dmom_v_adv_u_ext[nSpace],
6161 dmom_v_adv_v_ext[nSpace],
6162 dmom_v_adv_w_ext[nSpace],
6163 mom_w_adv_ext[nSpace],
6164 dmom_w_adv_u_ext[nSpace],
6165 dmom_w_adv_v_ext[nSpace],
6166 dmom_w_adv_w_ext[nSpace],
6167 mom_uu_diff_ten_ext[nSpace],
6168 mom_vv_diff_ten_ext[nSpace],
6169 mom_ww_diff_ten_ext[nSpace],
6170 mom_uv_diff_ten_ext[1],
6171 mom_uw_diff_ten_ext[1],
6172 mom_vu_diff_ten_ext[1],
6173 mom_vw_diff_ten_ext[1],
6174 mom_wu_diff_ten_ext[1],
6175 mom_wv_diff_ten_ext[1],
6176 mom_u_source_ext=0.0,
6177 mom_v_source_ext=0.0,
6178 mom_w_source_ext=0.0,
6180 dmom_u_ham_grad_p_ext[nSpace],
6181 dmom_u_ham_grad_u_ext[nSpace],
6183 dmom_v_ham_grad_p_ext[nSpace],
6184 dmom_v_ham_grad_v_ext[nSpace],
6186 dmom_w_ham_grad_p_ext[nSpace],
6187 dmom_w_ham_grad_w_ext[nSpace],
6188 dmom_u_adv_p_ext[nSpace],
6189 dmom_v_adv_p_ext[nSpace],
6190 dmom_w_adv_p_ext[nSpace],
6191 dflux_mass_u_ext=0.0,
6192 dflux_mass_v_ext=0.0,
6193 dflux_mass_w_ext=0.0,
6194 dflux_mom_u_adv_p_ext=0.0,
6195 dflux_mom_u_adv_u_ext=0.0,
6196 dflux_mom_u_adv_v_ext=0.0,
6197 dflux_mom_u_adv_w_ext=0.0,
6198 dflux_mom_v_adv_p_ext=0.0,
6199 dflux_mom_v_adv_u_ext=0.0,
6200 dflux_mom_v_adv_v_ext=0.0,
6201 dflux_mom_v_adv_w_ext=0.0,
6202 dflux_mom_w_adv_p_ext=0.0,
6203 dflux_mom_w_adv_u_ext=0.0,
6204 dflux_mom_w_adv_v_ext=0.0,
6205 dflux_mom_w_adv_w_ext=0.0,
6210 bc_mom_u_acc_ext=0.0,
6211 bc_dmom_u_acc_u_ext=0.0,
6212 bc_mom_v_acc_ext=0.0,
6213 bc_dmom_v_acc_v_ext=0.0,
6214 bc_mom_w_acc_ext=0.0,
6215 bc_dmom_w_acc_w_ext=0.0,
6216 bc_mass_adv_ext[nSpace],
6217 bc_dmass_adv_u_ext[nSpace],
6218 bc_dmass_adv_v_ext[nSpace],
6219 bc_dmass_adv_w_ext[nSpace],
6220 bc_mom_u_adv_ext[nSpace],
6221 bc_dmom_u_adv_u_ext[nSpace],
6222 bc_dmom_u_adv_v_ext[nSpace],
6223 bc_dmom_u_adv_w_ext[nSpace],
6224 bc_mom_v_adv_ext[nSpace],
6225 bc_dmom_v_adv_u_ext[nSpace],
6226 bc_dmom_v_adv_v_ext[nSpace],
6227 bc_dmom_v_adv_w_ext[nSpace],
6228 bc_mom_w_adv_ext[nSpace],
6229 bc_dmom_w_adv_u_ext[nSpace],
6230 bc_dmom_w_adv_v_ext[nSpace],
6231 bc_dmom_w_adv_w_ext[nSpace],
6232 bc_mom_uu_diff_ten_ext[nSpace],
6233 bc_mom_vv_diff_ten_ext[nSpace],
6234 bc_mom_ww_diff_ten_ext[nSpace],
6235 bc_mom_uv_diff_ten_ext[1],
6236 bc_mom_uw_diff_ten_ext[1],
6237 bc_mom_vu_diff_ten_ext[1],
6238 bc_mom_vw_diff_ten_ext[1],
6239 bc_mom_wu_diff_ten_ext[1],
6240 bc_mom_wv_diff_ten_ext[1],
6241 bc_mom_u_source_ext=0.0,
6242 bc_mom_v_source_ext=0.0,
6243 bc_mom_w_source_ext=0.0,
6244 bc_mom_u_ham_ext=0.0,
6245 bc_dmom_u_ham_grad_p_ext[nSpace],
6246 bc_dmom_u_ham_grad_u_ext[nSpace],
6247 bc_mom_v_ham_ext=0.0,
6248 bc_dmom_v_ham_grad_p_ext[nSpace],
6249 bc_dmom_v_ham_grad_v_ext[nSpace],
6250 bc_mom_w_ham_ext=0.0,
6251 bc_dmom_w_ham_grad_p_ext[nSpace],
6252 bc_dmom_w_ham_grad_w_ext[nSpace],
6253 fluxJacobian_p_p[nDOF_trial_element],
6254 fluxJacobian_p_u[nDOF_trial_element],
6255 fluxJacobian_p_v[nDOF_trial_element],
6256 fluxJacobian_p_w[nDOF_trial_element],
6257 fluxJacobian_u_p[nDOF_trial_element],
6258 fluxJacobian_u_u[nDOF_trial_element],
6259 fluxJacobian_u_v[nDOF_trial_element],
6260 fluxJacobian_u_w[nDOF_trial_element],
6261 fluxJacobian_v_p[nDOF_trial_element],
6262 fluxJacobian_v_u[nDOF_trial_element],
6263 fluxJacobian_v_v[nDOF_trial_element],
6264 fluxJacobian_v_w[nDOF_trial_element],
6265 fluxJacobian_w_p[nDOF_trial_element],
6266 fluxJacobian_w_u[nDOF_trial_element],
6267 fluxJacobian_w_v[nDOF_trial_element],
6268 fluxJacobian_w_w[nDOF_trial_element],
6269 jac_ext[nSpace*nSpace],
6271 jacInv_ext[nSpace*nSpace],
6272 boundaryJac[nSpace*(nSpace-1)],
6273 metricTensor[(nSpace-1)*(nSpace-1)],
6274 metricTensorDetSqrt,
6275 p_grad_trial_trace[nDOF_trial_element*nSpace],
6276 vel_grad_trial_trace[nDOF_trial_element*nSpace],
6278 p_test_dS[nDOF_test_element],
6279 vel_test_dS[nDOF_test_element],
6281 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
6282 vel_grad_test_dS[nDOF_trial_element*nSpace],
6286 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
6287 ck.calculateMapping_elementBoundary(eN,
6293 mesh_trial_trace_ref.data(),
6294 mesh_grad_trial_trace_ref.data(),
6295 boundaryJac_ref.data(),
6301 metricTensorDetSqrt,
6305 ck.calculateMappingVelocity_elementBoundary(eN,
6309 mesh_velocity_dof.data(),
6311 mesh_trial_trace_ref.data(),
6312 xt_ext,yt_ext,zt_ext,
6318 dS = metricTensorDetSqrt*dS_ref[kb];
6319 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
6320 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
6322 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
6323 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
6324 const double particle_eps = particle_epsFact * (useMetrics * h_phi + (1.0 - useMetrics) * elementDiameter[eN]);
6329 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
6332 p_ext = ebqe_p[ebNE_kb];
6333 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
6334 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
6337 for (
int I=0;I<nSpace;I++)
6338 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace+I];
6339 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
6340 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
6343 for (
int j=0;j<nDOF_trial_element;j++)
6346 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
6347 for (
int I=0;I<nSpace;I++)
6348 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
6353 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
6355 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
6356 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
6359 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
6363 double distance_to_omega_solid = 1e10;
6364 if (use_ball_as_particle == 1)
6366 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
6370 distance_to_omega_solid = ebq_global_phi_solid[ebN*nQuadraturePoints_elementBoundary+kb];
6372 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.),rhoSave, nuSave;
6381 elementDiameter[eN],
6382 smagorinskyConstant,
6383 turbulenceClosureModel,
6386 ebqe_vf_ext[ebNE_kb],
6387 ebqe_phi_ext[ebNE_kb],
6388 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
6389 distance_to_omega_solid,
6390 ebqe_kappa_phi_ext[ebNE_kb],
6402 ebqe_velocity_star[ebNE_kb_nSpace+0],
6403 ebqe_velocity_star[ebNE_kb_nSpace+1],
6404 ebqe_velocity_star[ebNE_kb_nSpace+1],
6428 mom_uu_diff_ten_ext,
6429 mom_vv_diff_ten_ext,
6430 mom_ww_diff_ten_ext,
6431 mom_uv_diff_ten_ext,
6432 mom_uw_diff_ten_ext,
6433 mom_vu_diff_ten_ext,
6434 mom_vw_diff_ten_ext,
6435 mom_wu_diff_ten_ext,
6436 mom_wv_diff_ten_ext,
6441 dmom_u_ham_grad_p_ext,
6442 dmom_u_ham_grad_u_ext,
6444 dmom_v_ham_grad_p_ext,
6445 dmom_v_ham_grad_v_ext,
6447 dmom_w_ham_grad_p_ext,
6448 dmom_w_ham_grad_w_ext,
6456 MATERIAL_PARAMETERS_AS_FUNCTION,
6457 ebqe_density_as_function[ebNE_kb],
6458 ebqe_dynamic_viscosity_as_function[ebNE_kb],
6461 use_ball_as_particle,
6464 ball_velocity.data(),
6465 ball_angular_velocity.data(),
6466 INT_BY_PARTS_PRESSURE);
6475 elementDiameter[eN],
6476 smagorinskyConstant,
6477 turbulenceClosureModel,
6480 bc_ebqe_vf_ext[ebNE_kb],
6481 bc_ebqe_phi_ext[ebNE_kb],
6482 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
6483 distance_to_omega_solid,
6484 ebqe_kappa_phi_ext[ebNE_kb],
6496 ebqe_velocity_star[ebNE_kb_nSpace+0],
6497 ebqe_velocity_star[ebNE_kb_nSpace+1],
6498 ebqe_velocity_star[ebNE_kb_nSpace+1],
6499 bc_eddy_viscosity_ext,
6501 bc_dmom_u_acc_u_ext,
6503 bc_dmom_v_acc_v_ext,
6505 bc_dmom_w_acc_w_ext,
6511 bc_dmom_u_adv_u_ext,
6512 bc_dmom_u_adv_v_ext,
6513 bc_dmom_u_adv_w_ext,
6515 bc_dmom_v_adv_u_ext,
6516 bc_dmom_v_adv_v_ext,
6517 bc_dmom_v_adv_w_ext,
6519 bc_dmom_w_adv_u_ext,
6520 bc_dmom_w_adv_v_ext,
6521 bc_dmom_w_adv_w_ext,
6522 bc_mom_uu_diff_ten_ext,
6523 bc_mom_vv_diff_ten_ext,
6524 bc_mom_ww_diff_ten_ext,
6525 bc_mom_uv_diff_ten_ext,
6526 bc_mom_uw_diff_ten_ext,
6527 bc_mom_vu_diff_ten_ext,
6528 bc_mom_vw_diff_ten_ext,
6529 bc_mom_wu_diff_ten_ext,
6530 bc_mom_wv_diff_ten_ext,
6531 bc_mom_u_source_ext,
6532 bc_mom_v_source_ext,
6533 bc_mom_w_source_ext,
6535 bc_dmom_u_ham_grad_p_ext,
6536 bc_dmom_u_ham_grad_u_ext,
6538 bc_dmom_v_ham_grad_p_ext,
6539 bc_dmom_v_ham_grad_v_ext,
6541 bc_dmom_w_ham_grad_p_ext,
6542 bc_dmom_w_ham_grad_w_ext,
6550 MATERIAL_PARAMETERS_AS_FUNCTION,
6551 ebqe_density_as_function[ebNE_kb],
6552 ebqe_dynamic_viscosity_as_function[ebNE_kb],
6555 use_ball_as_particle,
6558 ball_velocity.data(),
6559 ball_angular_velocity.data(),
6560 INT_BY_PARTS_PRESSURE);
6562 if (turbulenceClosureModel >= 3)
6564 const double turb_var_grad_0_dummy[2] = {0.,0.};
6565 const double c_mu = 0.09;
6574 ebqe_vf_ext[ebNE_kb],
6575 ebqe_phi_ext[ebNE_kb],
6578 ebqe_turb_var_0[ebNE_kb],
6579 ebqe_turb_var_1[ebNE_kb],
6580 turb_var_grad_0_dummy,
6582 mom_uu_diff_ten_ext,
6583 mom_vv_diff_ten_ext,
6584 mom_ww_diff_ten_ext,
6585 mom_uv_diff_ten_ext,
6586 mom_uw_diff_ten_ext,
6587 mom_vu_diff_ten_ext,
6588 mom_vw_diff_ten_ext,
6589 mom_wu_diff_ten_ext,
6590 mom_wv_diff_ten_ext,
6603 ebqe_vf_ext[ebNE_kb],
6604 ebqe_phi_ext[ebNE_kb],
6607 ebqe_turb_var_0[ebNE_kb],
6608 ebqe_turb_var_1[ebNE_kb],
6609 turb_var_grad_0_dummy,
6610 bc_eddy_viscosity_ext,
6611 bc_mom_uu_diff_ten_ext,
6612 bc_mom_vv_diff_ten_ext,
6613 bc_mom_ww_diff_ten_ext,
6614 bc_mom_uv_diff_ten_ext,
6615 bc_mom_uw_diff_ten_ext,
6616 bc_mom_vu_diff_ten_ext,
6617 bc_mom_vw_diff_ten_ext,
6618 bc_mom_wu_diff_ten_ext,
6619 bc_mom_wv_diff_ten_ext,
6620 bc_mom_u_source_ext,
6621 bc_mom_v_source_ext,
6622 bc_mom_w_source_ext);
6627 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
6628 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
6630 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
6631 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
6634 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
6635 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
6637 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
6638 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
6650 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
6651 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
6654 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
6655 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
6665 isDOFBoundary_u[ebNE_kb],
6666 isDOFBoundary_v[ebNE_kb],
6667 isDOFBoundary_w[ebNE_kb],
6668 isAdvectiveFluxBoundary_p[ebNE_kb],
6669 isAdvectiveFluxBoundary_u[ebNE_kb],
6670 isAdvectiveFluxBoundary_v[ebNE_kb],
6671 isAdvectiveFluxBoundary_w[ebNE_kb],
6672 dmom_u_ham_grad_p_ext[0],
6674 porosity_ext*dmom_u_acc_u_ext,
6683 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),
6684 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
6685 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
6686 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
6713 dflux_mom_u_adv_p_ext,
6714 dflux_mom_u_adv_u_ext,
6715 dflux_mom_u_adv_v_ext,
6716 dflux_mom_u_adv_w_ext,
6717 dflux_mom_v_adv_p_ext,
6718 dflux_mom_v_adv_u_ext,
6719 dflux_mom_v_adv_v_ext,
6720 dflux_mom_v_adv_w_ext,
6721 dflux_mom_w_adv_p_ext,
6722 dflux_mom_w_adv_u_ext,
6723 dflux_mom_w_adv_v_ext,
6724 dflux_mom_w_adv_w_ext,
6725 &ebqe_velocity_star[ebNE_kb_nSpace]);
6729 ck.calculateGScale(G,normal,h_penalty);
6730 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
6731 for (
int j=0;j<nDOF_trial_element;j++)
6733 register int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
6740 fluxJacobian_u_u[j] =
6741 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6743 ebqe_phi_ext[ebNE_kb],
6744 sdInfo_u_u_rowptr.data(),
6745 sdInfo_u_u_colind.data(),
6746 isDOFBoundary_u[ebNE_kb],
6747 isDiffusiveFluxBoundary_u[ebNE_kb],
6749 mom_uu_diff_ten_ext,
6750 vel_trial_trace_ref[ebN_local_kb_j],
6751 &vel_grad_trial_trace[j_nSpace],
6753 fluxJacobian_u_v[j]=
6754 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6756 ebqe_phi_ext[ebNE_kb],
6757 sdInfo_u_v_rowptr.data(),
6758 sdInfo_u_v_colind.data(),
6759 isDOFBoundary_v[ebNE_kb],
6760 isDiffusiveFluxBoundary_v[ebNE_kb],
6762 mom_uv_diff_ten_ext,
6763 vel_trial_trace_ref[ebN_local_kb_j],
6764 &vel_grad_trial_trace[j_nSpace],
6781 fluxJacobian_v_u[j]=
6782 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6784 ebqe_phi_ext[ebNE_kb],
6785 sdInfo_v_u_rowptr.data(),
6786 sdInfo_v_u_colind.data(),
6787 isDOFBoundary_u[ebNE_kb],
6788 isDiffusiveFluxBoundary_u[ebNE_kb],
6790 mom_vu_diff_ten_ext,
6791 vel_trial_trace_ref[ebN_local_kb_j],
6792 &vel_grad_trial_trace[j_nSpace],
6794 fluxJacobian_v_v[j]=
6795 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6797 ebqe_phi_ext[ebNE_kb],
6798 sdInfo_v_v_rowptr.data(),
6799 sdInfo_v_v_colind.data(),
6800 isDOFBoundary_v[ebNE_kb],
6801 isDiffusiveFluxBoundary_v[ebNE_kb],
6803 mom_vv_diff_ten_ext,
6804 vel_trial_trace_ref[ebN_local_kb_j],
6805 &vel_grad_trial_trace[j_nSpace],
6862 for (
int i=0;i<nDOF_test_element;i++)
6864 register int eN_i = eN*nDOF_test_element+i;
6865 for (
int j=0;j<nDOF_trial_element;j++)
6867 register int ebN_i_j = ebN*4*
nDOF_test_X_trial_element + i*nDOF_trial_element + j,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
6875 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*vel_test_dS[i]+
6876 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6877 isDiffusiveFluxBoundary_u[ebNE_kb],
6879 vel_trial_trace_ref[ebN_local_kb_j],
6881 sdInfo_u_u_rowptr.data(),
6882 sdInfo_u_u_colind.data(),
6883 mom_uu_diff_ten_ext,
6884 &vel_grad_test_dS[i*nSpace]);
6885 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] += fluxJacobian_u_v[j]*vel_test_dS[i]+
6886 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6887 isDiffusiveFluxBoundary_u[ebNE_kb],
6889 vel_trial_trace_ref[ebN_local_kb_j],
6891 sdInfo_u_v_rowptr.data(),
6892 sdInfo_u_v_colind.data(),
6893 mom_uv_diff_ten_ext,
6894 &vel_grad_test_dS[i*nSpace]);
6907 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] += fluxJacobian_v_u[j]*vel_test_dS[i]+
6908 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6909 isDiffusiveFluxBoundary_v[ebNE_kb],
6911 vel_trial_trace_ref[ebN_local_kb_j],
6913 sdInfo_v_u_rowptr.data(),
6914 sdInfo_v_u_colind.data(),
6915 mom_vu_diff_ten_ext,
6916 &vel_grad_test_dS[i*nSpace]);
6917 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] += fluxJacobian_v_v[j]*vel_test_dS[i]+
6918 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6919 isDiffusiveFluxBoundary_v[ebNE_kb],
6921 vel_trial_trace_ref[ebN_local_kb_j],
6923 sdInfo_v_v_rowptr.data(),
6924 sdInfo_v_v_colind.data(),
6925 mom_vv_diff_ten_ext,
6926 &vel_grad_test_dS[i*nSpace]);
6973 gf.useExact=useExact;
6974 gf_s.useExact=useExact;
6979 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
6980 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
6981 int nInteriorElementBoundaries_global = args.
scalar<
int>(
"nInteriorElementBoundaries_global");
6982 xt::pyarray<int>& interiorElementBoundariesArray = args.
array<
int>(
"interiorElementBoundariesArray");
6983 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
6984 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
6985 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
6986 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
6987 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
6988 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
6989 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
6990 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
6991 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
6992 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
6993 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
6994 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
6995 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
6996 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
6997 xt::pyarray<double>& vos_dof = args.
array<
double>(
"vos_dof");
6998 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
6999 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
7000 xt::pyarray<double>& velocityAverage = args.
array<
double>(
"velocityAverage");
7001 int permutations[nQuadraturePoints_elementBoundary];
7002 double xArray_left[nQuadraturePoints_elementBoundary*2],
7003 xArray_right[nQuadraturePoints_elementBoundary*2];
7004 for (
int i=0;i<nQuadraturePoints_elementBoundary;i++)
7006 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
7008 register int ebN = exteriorElementBoundariesArray[ebNE];
7009 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
7011 register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace,
7012 ebNE_kb_nSpace = ebNE*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace;
7013 velocityAverage[ebN_kb_nSpace+0]=ebqe_velocity[ebNE_kb_nSpace+0];
7014 velocityAverage[ebN_kb_nSpace+1]=ebqe_velocity[ebNE_kb_nSpace+1];
7017 for (
int ebNI = 0; ebNI < nInteriorElementBoundaries_global; ebNI++)
7019 register int ebN = interiorElementBoundariesArray[ebNI],
7020 left_eN_global = elementBoundaryElementsArray[ebN*2+0],
7021 left_ebN_element = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
7022 right_eN_global = elementBoundaryElementsArray[ebN*2+1],
7023 right_ebN_element = elementBoundaryLocalElementBoundariesArray[ebN*2+1],
7024 left_eN_nDOF_trial_element = left_eN_global*nDOF_trial_element,
7025 right_eN_nDOF_trial_element = right_eN_global*nDOF_trial_element;
7026 double jac[nSpace*nSpace],
7028 jacInv[nSpace*nSpace],
7029 boundaryJac[nSpace*(nSpace-1)],
7030 metricTensor[(nSpace-1)*(nSpace-1)],
7031 metricTensorDetSqrt,
7034 xt,yt,zt,integralScaling;
7036 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
7038 ck.calculateMapping_elementBoundary(left_eN_global,
7041 left_ebN_element*nQuadraturePoints_elementBoundary+kb,
7044 mesh_trial_trace_ref.data(),
7045 mesh_grad_trial_trace_ref.data(),
7046 boundaryJac_ref.data(),
7052 metricTensorDetSqrt,
7056 xArray_left[kb*2+0] = x;
7057 xArray_left[kb*2+1] = y;
7059 ck.calculateMapping_elementBoundary(right_eN_global,
7062 right_ebN_element*nQuadraturePoints_elementBoundary+kb,
7065 mesh_trial_trace_ref.data(),
7066 mesh_grad_trial_trace_ref.data(),
7067 boundaryJac_ref.data(),
7073 metricTensorDetSqrt,
7077 ck.calculateMappingVelocity_elementBoundary(left_eN_global,
7080 left_ebN_element*nQuadraturePoints_elementBoundary+kb,
7081 mesh_velocity_dof.data(),
7083 mesh_trial_trace_ref.data(),
7089 xArray_right[kb*2+0] = x;
7090 xArray_right[kb*2+1] = y;
7093 for (
int kb_left=0;kb_left<nQuadraturePoints_elementBoundary;kb_left++)
7095 double errorNormMin = 1.0;
7096 for (
int kb_right=0;kb_right<nQuadraturePoints_elementBoundary;kb_right++)
7098 double errorNorm=0.0;
7099 for (
int I=0;I<nSpace;I++)
7101 errorNorm += fabs(xArray_left[kb_left*2+I]
7103 xArray_right[kb_right*2+I]);
7105 if (errorNorm < errorNormMin)
7107 permutations[kb_right] = kb_left;
7108 errorNormMin = errorNorm;
7112 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
7114 register int ebN_kb_nSpace = ebN*nQuadraturePoints_elementBoundary*nSpace+kb*nSpace;
7115 register double u_left=0.0,
7125 register int left_kb = kb,
7126 right_kb = permutations[kb],
7127 left_ebN_element_kb_nDOF_test_element=(left_ebN_element*nQuadraturePoints_elementBoundary+left_kb)*nDOF_test_element,
7128 right_ebN_element_kb_nDOF_test_element=(right_ebN_element*nQuadraturePoints_elementBoundary+right_kb)*nDOF_test_element;
7132 ck.valFromDOF(vos_dof.data(),&vel_l2g[left_eN_nDOF_trial_element],&vel_trial_trace_ref[left_ebN_element_kb_nDOF_test_element],vos_left);
7133 ck.valFromDOF(u_dof.data(),&vel_l2g[left_eN_nDOF_trial_element],&vel_trial_trace_ref[left_ebN_element_kb_nDOF_test_element],u_left);
7134 ck.valFromDOF(v_dof.data(),&vel_l2g[left_eN_nDOF_trial_element],&vel_trial_trace_ref[left_ebN_element_kb_nDOF_test_element],v_left);
7137 ck.valFromDOF(vos_dof.data(),&vel_l2g[right_eN_nDOF_trial_element],&vel_trial_trace_ref[right_ebN_element_kb_nDOF_test_element],vos_right);
7138 ck.valFromDOF(u_dof.data(),&vel_l2g[right_eN_nDOF_trial_element],&vel_trial_trace_ref[right_ebN_element_kb_nDOF_test_element],u_right);
7139 ck.valFromDOF(v_dof.data(),&vel_l2g[right_eN_nDOF_trial_element],&vel_trial_trace_ref[right_ebN_element_kb_nDOF_test_element],v_right);
7144 velocityAverage[ebN_kb_nSpace+0]=0.5*(u_left + u_right);
7145 velocityAverage[ebN_kb_nSpace+1]=0.5*(v_left + v_right);
7154 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
7155 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
7156 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
7157 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
7158 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
7159 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
7160 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
7161 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
7162 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
7163 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
7164 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
7165 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
7166 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
7167 xt::pyarray<double>& isBoundary_1D = args.
array<
double>(
"isBoundary_1D");
7174 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
7176 register int ebN = exteriorElementBoundariesArray[ebNE],
7177 eN = elementBoundaryElementsArray[ebN*2+0],
7178 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
7179 eN_nDOF_trial_element = eN*nDOF_trial_element;
7181 elementIsBoundary[nDOF_test_element];
7182 const double* elementResidual_w(NULL);
7183 for (
int i=0;i<nDOF_test_element;i++)
7184 elementIsBoundary[i]=0.0;
7185 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
7187 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
7188 ebNE_kb_nSpace = ebNE_kb*nSpace,
7189 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
7190 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
7192 jac_ext[nSpace*nSpace],
7194 jacInv_ext[nSpace*nSpace],
7195 boundaryJac[nSpace*(nSpace-1)],
7196 metricTensor[(nSpace-1)*(nSpace-1)],
7197 metricTensorDetSqrt,
7198 dS, vel_test_dS[nDOF_test_element],
7199 normal[2],x_ext,y_ext,z_ext;
7201 ck.calculateMapping_elementBoundary(eN,
7207 mesh_trial_trace_ref.data(),
7208 mesh_grad_trial_trace_ref.data(),
7209 boundaryJac_ref.data(),
7215 metricTensorDetSqrt,
7219 dS = metricTensorDetSqrt*dS_ref[kb];
7221 for (
int j=0;j<nDOF_trial_element;j++)
7222 vel_test_dS[j] = fabs(vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j])*dS;
7226 for (
int i=0;i<nDOF_test_element;i++)
7227 elementIsBoundary[i] += vel_test_dS[i];
7232 for (
int i=0;i<nDOF_test_element;i++)
7234 int eN_i = eN*nDOF_test_element+i;
7235 isBoundary_1D[vel_l2g[eN_i]] += elementIsBoundary[i];
7242 int nQuadraturePoints_elementIn,
7243 int nDOF_mesh_trial_elementIn,
7244 int nDOF_trial_elementIn,
7245 int nDOF_test_elementIn,
7246 int nQuadraturePoints_elementBoundaryIn,
7251 double packFraction,
7262 double angFriction,
double vos_limiter,
double mu_fr_limiter )
7264 cppRANS3PF2D_base *rvalue = proteus::chooseAndAllocateDiscretization2D<cppRANS3PF2D_base, cppRANS3PF2D, CompKernel>(nSpaceIn,
7265 nQuadraturePoints_elementIn,
7266 nDOF_mesh_trial_elementIn,
7267 nDOF_trial_elementIn,
7268 nDOF_test_elementIn,
7269 nQuadraturePoints_elementBoundaryIn,
7285 angFriction, vos_limiter, mu_fr_limiter );