1692 double NONCONSERVATIVE_FORM = args.
scalar<
double>(
"NONCONSERVATIVE_FORM");
1693 double MOMENTUM_SGE = args.
scalar<
double>(
"MOMENTUM_SGE");
1694 double PRESSURE_SGE = args.
scalar<
double>(
"PRESSURE_SGE");
1695 double VELOCITY_SGE = args.
scalar<
double>(
"VELOCITY_SGE");
1696 double PRESSURE_PROJECTION_STABILIZATION = args.
scalar<
double>(
"PRESSURE_PROJECTION_STABILIZATION");
1697 xt::pyarray<double>& numerical_viscosity = args.
array<
double>(
"numerical_viscosity");
1698 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1699 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1700 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1701 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1702 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1703 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1704 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
1705 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1706 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1707 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1708 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1709 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1710 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1711 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1712 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1713 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1714 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1715 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1716 xt::pyarray<double>& xb_ref = args.
array<
double>(
"xb_ref");
1717 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1718 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1719 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1720 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1721 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1722 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1723 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1724 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1725 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1726 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1727 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1728 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1729 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1730 xt::pyarray<double>& elementBoundaryDiameter = args.
array<
double>(
"elementBoundaryDiameter");
1731 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1732 double hFactor = args.
scalar<
double>(
"hFactor");
1733 int nElements_global = args.
scalar<
int>(
"nElements_global");
1734 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1735 double useRBLES = args.
scalar<
double>(
"useRBLES");
1736 double useMetrics = args.
scalar<
double>(
"useMetrics");
1737 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1738 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1739 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1740 double sigma = args.
scalar<
double>(
"sigma");
1745 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1746 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1747 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1748 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1749 double C_dc = args.
scalar<
double>(
"C_dc");
1750 double C_b = args.
scalar<
double>(
"C_b");
1751 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1752 xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1753 const xt::pyarray<double>& eps_porous = args.
array<
double>(
"eps_porous");
1754 xt::pyarray<double>& phi_porous = args.
array<
double>(
"phi_porous");
1755 const xt::pyarray<double>& q_velocity_porous = args.
array<
double>(
"q_velocity_porous");
1756 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
1757 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1758 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1759 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1760 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1761 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1762 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1763 const double LAG_LES = args.
scalar<
double>(
"LAG_LES");
1764 xt::pyarray<double> & q_eddy_viscosity = args.
array<
double>(
"q_eddy_viscosity");
1765 xt::pyarray<double> & q_eddy_viscosity_last = args.
array<
double>(
"q_eddy_viscosity_last");
1766 xt::pyarray<double> & ebqe_eddy_viscosity = args.
array<
double>(
"ebqe_eddy_viscosity");
1767 xt::pyarray<double> & ebqe_eddy_viscosity_last = args.
array<
double>(
"ebqe_eddy_viscosity_last");
1768 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1769 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1770 xt::pyarray<int>& rp_l2g = args.
array<
int>(
"rp_l2g");
1771 xt::pyarray<int>& rvel_l2g = args.
array<
int>(
"rvel_l2g");
1772 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1773 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1774 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1775 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1776 xt::pyarray<double>& p_old_dof = args.
array<
double>(
"p_old_dof");
1777 xt::pyarray<double>& u_old_dof = args.
array<
double>(
"u_old_dof");
1778 xt::pyarray<double>& v_old_dof = args.
array<
double>(
"v_old_dof");
1779 xt::pyarray<double>& w_old_dof = args.
array<
double>(
"w_old_dof");
1780 xt::pyarray<double>& g = args.
array<
double>(
"g");
1781 const double useVF = args.
scalar<
double>(
"useVF");
1782 xt::pyarray<double>& q_rho = args.
array<
double>(
"q_rho");
1783 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1784 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1785 xt::pyarray<double>& phi_nodes = args.
array<
double>(
"phi_nodes");
1786 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1787 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1788 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1789 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1790 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1791 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1792 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1793 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1794 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1795 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
1796 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
1797 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1798 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1799 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1800 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1801 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1802 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1803 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1804 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1805 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1806 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1807 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1808 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1809 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1810 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1811 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1812 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1813 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1814 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1815 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1816 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1817 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1818 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1819 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1820 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1821 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1822 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1823 int offset_p = args.
scalar<
int>(
"offset_p");
1824 int offset_u = args.
scalar<
int>(
"offset_u");
1825 int offset_v = args.
scalar<
int>(
"offset_v");
1826 int offset_w = args.
scalar<
int>(
"offset_w");
1827 int stride_p = args.
scalar<
int>(
"stride_p");
1828 int stride_u = args.
scalar<
int>(
"stride_u");
1829 int stride_v = args.
scalar<
int>(
"stride_v");
1830 int stride_w = args.
scalar<
int>(
"stride_w");
1831 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1832 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1833 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1834 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
1835 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1836 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1837 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1838 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1839 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1840 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1841 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1842 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1843 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
1844 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1845 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1846 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1847 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1848 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1849 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1850 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1851 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1852 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1853 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1854 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1855 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1856 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1857 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1858 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1859 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1860 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1861 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1862 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1863 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1864 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1865 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1866 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1867 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1868 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1869 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1870 xt::pyarray<double>& q_u_0 = args.
array<
double>(
"q_u_0");
1871 xt::pyarray<double>& q_u_1 = args.
array<
double>(
"q_u_1");
1872 xt::pyarray<double>& q_u_2 = args.
array<
double>(
"q_u_2");
1873 xt::pyarray<double>& q_u_3 = args.
array<
double>(
"q_u_3");
1874 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1875 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1876 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1877 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1878 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
1879 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1880 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1881 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1882 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1883 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1884 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1885 xt::pyarray<double>& velocityError = args.
array<
double>(
"velocityError");
1886 xt::pyarray<double>& velocityErrorNodal = args.
array<
double>(
"velocityErrorNodal");
1887 xt::pyarray<double>& forcex = args.
array<
double>(
"forcex");
1888 xt::pyarray<double>& forcey = args.
array<
double>(
"forcey");
1889 xt::pyarray<double>& forcez = args.
array<
double>(
"forcez");
1890 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
1891 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
1892 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
1893 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
1894 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
1895 xt::pyarray<double>& ball_density = args.
array<
double>(
"ball_density");
1896 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
1897 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
1898 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
1899 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
1900 xt::pyarray<double>& ebqe_phi_s = args.
array<
double>(
"ebqe_phi_s");
1901 xt::pyarray<double>& ebq_global_grad_phi_s = args.
array<
double>(
"ebq_global_grad_phi_s");
1902 xt::pyarray<double>& ebq_particle_velocity_s = args.
array<
double>(
"ebq_particle_velocity_s");
1903 int nParticles = args.
scalar<
int>(
"nParticles");
1904 xt::pyarray<double>& particle_netForces = args.
array<
double>(
"particle_netForces");
1905 xt::pyarray<double>& particle_netMoments = args.
array<
double>(
"particle_netMoments");
1906 xt::pyarray<double>& particle_surfaceArea = args.
array<
double>(
"particle_surfaceArea");
1907 xt::pyarray<double>& particle_surfaceArea_projected = args.
array<
double>(
"particle_surfaceArea_projected");
1908 xt::pyarray<double>& projection_direction = args.
array<
double>(
"projection_direction");
1909 xt::pyarray<double>& particle_volume = args.
array<
double>(
"particle_volume");
1910 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1911 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
1912 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
1913 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
1914 double particle_beta = args.
scalar<
double>(
"particle_beta");
1915 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
1916 double ghost_penalty_constant = args.
scalar<
double>(
"ghost_penalty_constant");
1917 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
1918 xt::pyarray<double>& distance_to_solids = args.
array<
double>(
"distance_to_solids");
1919 bool useExact = args.
scalar<
int>(
"useExact");
1920 xt::pyarray<double>& isActiveR = args.
array<
double>(
"isActiveR");
1921 xt::pyarray<double>& isActiveDOF_p = args.
array<
double>(
"isActiveDOF_p");
1922 xt::pyarray<double>& isActiveDOF_vel = args.
array<
double>(
"isActiveDOF_vel");
1923 const bool normalize_pressure = args.
scalar<
int>(
"normalize_pressure");
1924 xt::pyarray<double>& errors = args.
array<
double>(
"errors");
1925 xt::pyarray<double>& ball_u = args.
array<
double>(
"ball_u");
1926 xt::pyarray<double>& ball_v = args.
array<
double>(
"ball_v");
1927 xt::pyarray<int>& isActiveElement = args.
array<
int>(
"isActiveElement");
1928 xt::pyarray<int>& isActiveElement_last = args.
array<
int>(
"isActiveElement_last");
1929 logEvent(
"Entered mprans calculateResidual",6);
1937 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
1941 double p_dv=0.0,pa_dv=0.0,total_volume=0.0,total_surface_area=0.0,total_flux=0.0;
1942 double mesh_volume_conservation=0.0,
1943 mesh_volume_conservation_weak=0.0,
1944 mesh_volume_conservation_err_max=0.0,
1945 mesh_volume_conservation_err_max_weak=0.0,
1947 &p_L1=errors(0,0),&u_L1=errors(0,1),&v_L1=errors(0,2),&w_L1=errors(0,2),&velocity_L1=errors(0,4),
1948 &p_L2=errors(1,0),&u_L2=errors(1,1),&v_L2=errors(1,2),&w_L2=errors(1,2),&velocity_L2=errors(1,4),
1949 &p_LI=errors(2,0),&u_LI=errors(2,1),&v_LI=errors(2,2),&w_LI=errors(2,2),&velocity_LI=errors(2,4);
1950 p_L1=0.0; u_L1=0.0; v_L1=0.0; w_L1=0.0; velocity_L1=0.0;
1951 p_L2=0.0; u_L2=0.0; v_L2=0.0; w_L2=0.0; velocity_L2=0.0;
1952 p_LI=0.0; u_LI=0.0; v_LI=0.0; w_LI=0.0; velocity_LI=0.0;
1953 double globalConservationError=0.0;
1958 for(
int eN=0;eN<nElements_global;eN++)
1961 double elementResidual_p[nDOF_test_element],elementResidual_p_check[nDOF_test_element],elementResidual_mesh[nDOF_test_element],
1962 elementResidual_u[nDOF_v_test_element],
1963 elementResidual_v[nDOF_v_test_element],
1964 pelementResidual_u[nDOF_v_test_element],
1965 pelementResidual_v[nDOF_v_test_element],
1966 velocityErrorElement[nDOF_v_test_element],
1968 bool element_active=
false;
1969 isActiveElement[eN]=0;
1970 const double* elementResidual_w(NULL);
1971 double mesh_volume_conservation_element=0.0,
1972 mesh_volume_conservation_element_weak=0.0;
1973 int particle_index=0;
1974 for (
int i=0;i<nDOF_test_element;i++)
1976 int eN_i = eN*nDOF_test_element+i;
1977 elementResidual_p_save.data()[eN_i]=0.0;
1978 elementResidual_mesh[i]=0.0;
1979 elementResidual_p[i]=0.0;
1980 elementResidual_p_check[i]=0.0;
1982 for (
int i=0;i<nDOF_v_test_element;i++)
1984 elementResidual_u[i]=0.0;
1985 elementResidual_v[i]=0.0;
1986 pelementResidual_u[i]=0.0;
1987 pelementResidual_v[i]=0.0;
1988 velocityErrorElement[i]=0.0;
1991 if(use_ball_as_particle==1 && nParticles > 0)
1993 double min_d = 1e10;
1994 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1997 mesh_dof.data()[3*mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]+0],
1998 mesh_dof.data()[3*mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]+1],
1999 mesh_dof.data()[3*mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]+2],
2000 phi_solid_nodes.data()[mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]]);
2001 if (phi_solid_nodes.data()[mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]] < min_d)
2003 min_d = phi_solid_nodes.data()[mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]];
2004 particle_index = index;
2012 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
2013 for (
int j=0;j<nDOF_mesh_trial_element;j++)
2015 int eN_j = eN*nDOF_mesh_trial_element+j;
2016 element_phi[j] = phi_nodes.data()[p_l2g.data()[eN_j]];
2017 element_phi_s[j] = phi_solid_nodes.data()[p_l2g.data()[eN_j]];
2019 double element_nodes[nDOF_mesh_trial_element*3];
2020 for (
int i=0;i<nDOF_mesh_trial_element;i++)
2022 int eN_i=eN*nDOF_mesh_trial_element+i;
2023 for(
int I=0;I<3;I++)
2024 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
2026 int icase_s =
gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
2029 element_active=
true;
2030 isActiveElement[eN]=1;
2032 for (
int ebN_element=0;ebN_element < nDOF_mesh_trial_element; ebN_element++)
2034 const int ebN = elementBoundariesArray.data()[eN*nDOF_mesh_trial_element+ebN_element];
2037 if (elementBoundaryElementsArray[ebN*2+1] != -1 && element_phi_s[(ebN_element+1)%nDOF_mesh_trial_element]*element_phi_s[(ebN_element+2)%nDOF_mesh_trial_element] <= 0.0)
2040 if (elementBoundaryElementsArray[ebN*2 + 0] == eN)
2045 else if (icase_s == 1)
2047 element_active=
true;
2048 isActiveElement[eN]=1;
2051 int icase_p =
gf_p.
calculate(element_phi, element_nodes, x_ref.data(), -
rho_1*g.data()[1], -
rho_0*g.data()[1],
false,
true);
2054 int icase_p =
gf_p.
calculate(element_phi, element_nodes, x_ref.data(), 1.,1.,
false,
false);
2055 int icase =
gf.
calculate(element_phi, element_nodes, x_ref.data(), 1.,1.,
false,
false);
2060 for (
int ebN_element=0;ebN_element < nDOF_mesh_trial_element; ebN_element++)
2062 const int ebN = elementBoundariesArray.data()[eN*nDOF_mesh_trial_element+ebN_element];
2071 double numDiffMax=0.0;
2072 for(
int fluid_phase=0;fluid_phase < 2 - abs(icase);fluid_phase++)
2074 for(
int k=0;k<nQuadraturePoints_element;k++)
2077 int eN_k = eN*nQuadraturePoints_element+k,
2078 eN_k_nSpace = eN_k*nSpace,
2080 eN_nDOF_trial_element = eN*nDOF_trial_element,
2081 eN_nDOF_v_trial_element = eN*nDOF_v_trial_element;
2082 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
2084 p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0,
2112 mom_uu_diff_ten[nSpace]=
ZEROVEC,
2113 mom_vv_diff_ten[nSpace]=
ZEROVEC,
2114 mom_ww_diff_ten[nSpace]=
ZEROVEC,
2125 dmom_u_ham_grad_p[nSpace]=
ZEROVEC,
2126 dmom_u_ham_grad_u[nSpace]=
ZEROVEC,
2127 dmom_u_ham_grad_v[nSpace]=
ZEROVEC,
2132 dmom_v_ham_grad_p[nSpace]=
ZEROVEC,
2133 dmom_v_ham_grad_u[nSpace]=
ZEROVEC,
2134 dmom_v_ham_grad_v[nSpace]=
ZEROVEC,
2139 dmom_w_ham_grad_p[nSpace]=
ZEROVEC,
2140 dmom_w_ham_grad_w[nSpace]=
ZEROVEC,
2154 Lstar_u_p[nDOF_test_element],
2155 Lstar_v_p[nDOF_test_element],
2156 Lstar_w_p[nDOF_test_element],
2157 Lstar_u_u[nDOF_v_test_element],
2158 Lstar_v_v[nDOF_v_test_element],
2159 Lstar_w_w[nDOF_v_test_element],
2160 Lstar_p_u[nDOF_v_test_element],
2161 Lstar_p_v[nDOF_v_test_element],
2162 Lstar_p_w[nDOF_v_test_element],
2167 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
2168 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2171 jacInv[nSpace*nSpace],
2172 p_trial[nDOF_trial_element], vel_trial[nDOF_v_trial_element],
2173 p_grad_trial_ib[nDOF_trial_element*nSpace], vel_grad_trial_ib[nDOF_v_trial_element*nSpace],
2174 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_v_trial_element*nSpace],
2175 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_v_test_element],
2176 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_v_test_element*nSpace],
2183 dmom_u_source[nSpace]=
ZEROVEC,
2184 dmom_v_source[nSpace]=
ZEROVEC,
2185 dmom_w_source[nSpace]=
ZEROVEC,
2187 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace]=
ZEROVEC,dmom_adv_sge[nSpace]=
ZEROVEC,dmom_ham_grad_sge[nSpace]=
ZEROVEC,
2193 dmom_u_source_s[nSpace]=
ZEROVEC,
2194 dmom_v_source_s[nSpace]=
ZEROVEC,
2195 dmom_w_source_s[nSpace]=
ZEROVEC,
2199 dmom_u_adv_u_s[nSpace]=
ZEROVEC,
2200 dmom_v_adv_v_s[nSpace]=
ZEROVEC,
2201 dmom_w_adv_w_s[nSpace]=
ZEROVEC,
2203 dmom_u_ham_grad_u_s[nSpace]=
ZEROVEC,
2204 dmom_u_ham_grad_v_s[nSpace]=
ZEROVEC,
2209 dmom_v_ham_grad_u_s[nSpace]=
ZEROVEC,
2210 dmom_v_ham_grad_v_s[nSpace]=
ZEROVEC,
2215 dmom_w_ham_grad_w_s[nSpace]=
ZEROVEC,
2227 ck.calculateMapping_element(eN,
2231 mesh_trial_ref.data(),
2232 mesh_grad_trial_ref.data(),
2237 ck.calculateH_element(eN,
2239 nodeDiametersArray.data(),
2241 mesh_trial_ref.data(),
2244 ck.calculateMappingVelocity_element(eN,
2246 mesh_velocity_dof.data(),
2248 mesh_trial_ref.data(),
2253 dV = fabs(jacDet)*dV_ref.data()[k];
2254 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2257 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2258 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2261 ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
2262 ck_v.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_v_trial_element*nSpace],jacInv,vel_grad_trial);
2263 for (
int i=0; i < nDOF_trial_element; i++)
2265 p_trial[i] = p_trial_ref.data()[k*nDOF_trial_element + i];
2266 p_grad_trial_ib[i*nSpace + 0] = p_grad_trial[i*nSpace+0];
2267 p_grad_trial_ib[i*nSpace + 1] = p_grad_trial[i*nSpace+1];
2269 for (
int i=0; i < nDOF_v_trial_element; i++)
2271 vel_trial[i] = vel_trial_ref.data()[k*nDOF_v_trial_element + i];
2272 vel_grad_trial_ib[i*nSpace + 0] = vel_grad_trial[i*nSpace+0];
2273 vel_grad_trial_ib[i*nSpace + 1] = vel_grad_trial[i*nSpace+1];
2278 for (
int i=0; i < nDOF_trial_element; i++)
2280 if (fluid_phase == 0)
2282 if (not std::isnan(
gf_p.
VA(i)))
2285 p_grad_trial_ib[i*nSpace + 0] =
gf_p.
VA_x(i);
2286 p_grad_trial_ib[i*nSpace + 1] =
gf_p.
VA_y(i);
2291 if (not std::isnan(
gf_p.
VB(i)))
2294 p_grad_trial_ib[i*nSpace + 0] =
gf_p.
VB_x(i);
2295 p_grad_trial_ib[i*nSpace + 1] =
gf_p.
VB_y(i);
2299 if(nDOF_v_trial_element == nDOF_trial_element)
2301 for (
int vi=0; vi < nDOF_v_trial_element; vi++)
2303 if (fluid_phase == 0)
2305 if (not std::isnan(
gf.
VA(vi)))
2307 vel_trial[vi] =
gf.
VA(vi);
2308 vel_grad_trial_ib[vi*nSpace + 0] =
gf.
VA_x(vi);
2309 vel_grad_trial_ib[vi*nSpace + 1] =
gf.
VA_y(vi);
2314 if (not std::isnan(
gf.
VB(vi)))
2316 vel_trial[vi] =
gf.
VB(vi);
2317 vel_grad_trial_ib[vi*nSpace + 0] =
gf.
VB_x(vi);
2318 vel_grad_trial_ib[vi*nSpace + 1] =
gf.
VB_y(vi);
2326 for (
int vi=0; vi < nDOF_v_trial_element; vi++)
2329 if (fabs(p_trial_ref.data()[k*nDOF_trial_element + vi] - p_trial[vi]) > 1.0e-8)
2331 for (
int vj=0; vj < nDOF_trial_element; vj++)
2332 std::cout<<
"Trial "<<p_trial_ref.data()[k*nDOF_trial_element + vj]<<
'\t'<<
gf_p.
VA(vj)<<
'\t'<<
gf_p.
VB(vj)<<std::endl;
2335 if (fabs(p_grad_trial[vi*nSpace + 0] - p_grad_trial_ib[vi*nSpace+0]) > 1.0e-8)
2337 for (
int vj=0; vj < nDOF_trial_element; vj++)
2338 std::cout<<
"Grad Trial x"<<p_grad_trial[vj*nSpace + 0]<<
'\t'<<
gf_p.
VA_x(vj)<<
'\t'<<
gf_p.
VB_x(vj)<<std::endl;
2341 if (fabs(p_grad_trial[vi*nSpace + 1] - p_grad_trial_ib[vi*nSpace+1]) > 1.0e-8)
2343 for (
int vj=0; vj < nDOF_trial_element; vj++)
2344 std::cout<<
"Grad Trial y "<<p_grad_trial[vj*nSpace + 1]<<
'\t'<<
gf_p.
VA_y(vj)<<
'\t'<<
gf_p.
VB_y(vj)<<std::endl;
2348 if (fabs(vel_trial_ref.data()[k*nDOF_v_trial_element + vi] - vel_trial[vi]) > 1.0e-8)
2350 for (
int vj=0; vj < nDOF_v_trial_element; vj++)
2351 std::cout<<
"Trial "<<vel_trial_ref.data()[k*nDOF_v_trial_element + vj]<<
'\t'<<
gf.
VA(vj)<<
'\t'<<
gf.
VB(vj)<<std::endl;
2354 if (fabs(vel_grad_trial[vi*nSpace + 0] - vel_grad_trial_ib[vi*nSpace+0]) > 1.0e-8)
2356 for (
int vj=0; vj < nDOF_v_trial_element; vj++)
2357 std::cout<<
"Grad Trial x"<<vel_grad_trial[vj*nSpace + 0]<<
'\t'<<
gf.
VA_x(vj)<<
'\t'<<
gf.
VB_x(vj)<<std::endl;
2360 if (fabs(vel_grad_trial[vi*nSpace + 1] - vel_grad_trial_ib[vi*nSpace+1]) > 1.0e-8)
2362 for (
int vj=0; vj < nDOF_v_trial_element; vj++)
2363 std::cout<<
"Grad Trial y "<<vel_grad_trial[vj*nSpace + 1]<<
'\t'<<
gf.
VA_y(vj)<<
'\t'<<
gf.
VB_y(vj)<<std::endl;
2373 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_trial,p);
2374 ck_v.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_trial,
u);
2375 ck_v.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_trial,
v);
2376 ck.valFromDOF(p_old_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_trial,p_old);
2377 ck_v.valFromDOF(u_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_trial,u_old);
2378 ck_v.valFromDOF(v_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_trial,v_old);
2380 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_ib,grad_p);
2381 ck_v.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_u);
2382 ck_v.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_v);
2383 ck.gradFromDOF(p_old_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_ib,grad_p_old);
2384 ck_v.gradFromDOF(u_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_u_old);
2385 ck_v.gradFromDOF(v_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_v_old);
2387 if (PRESSURE_PROJECTION_STABILIZATION)
2388 ck.DOFaverage(p_dof.data(), &p_l2g.data()[eN_nDOF_trial_element],p_element_avg);
2391 for (
int j=0;j<nDOF_test_element;j++)
2393 p_test_dV[j] = p_trial[j]*dV;
2394 for (
int I=0;I<nSpace;I++)
2396 p_grad_test_dV[j*nSpace+I] = p_grad_trial_ib[j*nSpace+I]*dV;
2400 for (
int j=0;j<nDOF_v_test_element;j++)
2402 vel_test_dV[j] = vel_trial[j]*dV;
2403 for (
int I=0;I<nSpace;I++)
2405 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial_ib[j*nSpace+I]*dV;
2409 for (
int j=0;j<nDOF_test_element;j++)
2411 p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
2412 for (
int I=0;I<nSpace;I++)
2414 p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;
2418 for (
int j=0;j<nDOF_v_test_element;j++)
2420 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_v_trial_element+j]*dV;
2421 for (
int I=0;I<nSpace;I++)
2423 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
2428 double div_mesh_velocity=0.0;
2429 for (
int j=0;j<nDOF_trial_element;j++)
2431 int eN_j=eN*nDOF_trial_element+j;
2432 div_mesh_velocity +=
2433 mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+0]*p_grad_trial[j*nSpace+0] +
2434 mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+1]*p_grad_trial[j*nSpace+1];
2436 mesh_volume_conservation_element += (alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity)*dV;
2437 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last.data()[eN_k])/dV;
2439 porosity = q_porosity.data()[eN_k];
2441 q_velocity.data()[eN_k_nSpace+0]=
u;
2442 q_velocity.data()[eN_k_nSpace+1]=
v;
2443 q_x.data()[eN_k_3d + 0] = x;
2444 q_x.data()[eN_k_3d + 1] = y;
2445 double ball_n[nSpace];
2446 if (use_ball_as_particle == 1 && nParticles > 0)
2448 int ball_index=
get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x,y,
z,distance_to_solids.data()[eN_k]);
2449 get_normal_to_ith_ball(nParticles, ball_center.data(), ball_radius.data(),ball_index,x,y,
z,ball_n[0],ball_n[1]);
2456 phi_solid.data()[eN_k] = distance_to_solids.data()[eN_k];
2457 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2461 const double H_s =
gf_s.
H(particle_eps,phi_solid.data()[eN_k]);
2462 const double D_s =
gf_s.
D(particle_eps,phi_solid.data()[eN_k]);
2464 double p_e = q_u_0.data()[eN_k] - p,
2465 u_e = q_u_1.data()[eN_k] -
u,
2466 v_e = q_u_2.data()[eN_k] -
v,
2467 velocity_e=sqrt(u_e*u_e + v_e*v_e);
2480 if (fluid_phase == 0)
2491 else if (icase == -1)
2496 else if (icase == 1)
2506 double H = (1.0-useVF)*
gf.
H(eps_rho,
phi[eN_k]) + useVF*fmin(1.0,fmax(0.0,vf[eN_k]));
2507 double ImH = (1.0-useVF)*
gf.
ImH(eps_rho,
phi[eN_k]) + useVF*(1.0-fmin(1.0,fmax(0.0,vf[eN_k])));
2516 elementDiameter.data()[eN],
2517 smagorinskyConstant,
2518 turbulenceClosureModel,
2523 &normal_phi.data()[eN_k_nSpace],
2524 kappa_phi.data()[eN_k],
2527 phi_solid.data()[eN_k],
2546 q_eddy_viscosity.data()[eN_k],
2547 q_eddy_viscosity_last.data()[eN_k],
2600 forcex.data()[eN_k],
2601 forcey.data()[eN_k],
2602 forcez.data()[eN_k]);
2603 q_rho.data()[eN_k] = rho;
2605 mass_source = q_mass_source.data()[eN_k];
2612 q_dragAlpha.data()[eN_k],
2613 q_dragBeta.data()[eN_k],
2626 q_velocity_sge.data()[eN_k_nSpace+0],
2627 q_velocity_sge.data()[eN_k_nSpace+1],
2628 q_velocity_sge.data()[eN_k_nSpace+1],
2629 eps_porous.data()[elementFlags.data()[eN]],
2630 phi_porous.data()[eN_k],
2631 q_velocity_porous.data()[eN_k_nSpace+0],
2632 q_velocity_porous.data()[eN_k_nSpace+1],
2633 q_velocity_porous.data()[eN_k_nSpace+1],
2641 if (turbulenceClosureModel >= 3)
2643 const double c_mu = 0.09;
2645 turbulenceClosureModel,
2657 q_turb_var_0.data()[eN_k],
2658 q_turb_var_1.data()[eN_k],
2659 &q_turb_var_grad_0.data()[eN_k_nSpace],
2660 q_eddy_viscosity.data()[eN_k],
2677 if (NONCONSERVATIVE_FORM > 0.0)
2679 mom_u_ham -= MOVING_DOMAIN*dmom_u_acc_u*(grad_u[0]*
xt + grad_u[1]*yt);
2680 dmom_u_ham_grad_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
2681 dmom_u_ham_grad_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
2685 mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*
xt;
2686 mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt;
2687 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
2688 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
2691 if (NONCONSERVATIVE_FORM > 0.0)
2693 mom_v_ham -= MOVING_DOMAIN*dmom_v_acc_v*(grad_v[0]*
xt + grad_v[1]*yt);
2694 dmom_v_ham_grad_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
2695 dmom_v_ham_grad_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
2699 mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*
xt;
2700 mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt;
2701 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
2702 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
2708 if (q_dV_last.data()[eN_k] <= -100)
2709 q_dV_last.data()[eN_k] = dV;
2710 q_dV.data()[eN_k] = dV;
2712 q_mom_u_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
2718 q_mom_v_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
2724 if (NONCONSERVATIVE_FORM > 0.0)
2726 mom_u_acc_t *= dmom_u_acc_u;
2727 mom_v_acc_t *= dmom_v_acc_v;
2733 pdeResidual_p =
ck.Advection_strong(dmass_adv_u,grad_u) +
2734 ck.Advection_strong(dmass_adv_v,grad_v) +
2735 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity) +
2736 ck.Reaction_strong(mass_source);
2738 if (NONCONSERVATIVE_FORM > 0.0)
2740 dmom_adv_sge[0] = 0.0;
2741 dmom_adv_sge[1] = 0.0;
2742 dmom_ham_grad_sge[0] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
2743 dmom_ham_grad_sge[1] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
2747 dmom_adv_sge[0] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
2748 dmom_adv_sge[1] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
2749 dmom_ham_grad_sge[0] = 0.0;
2750 dmom_ham_grad_sge[1] = 0.0;
2752 double mv_tau[nSpace]=
ZEROVEC;
2753 mv_tau[0] = dmom_adv_sge[0] + dmom_ham_grad_sge[0];
2754 mv_tau[1] = dmom_adv_sge[1] + dmom_ham_grad_sge[1];
2756 pdeResidual_u =
ck.Mass_strong(mom_u_acc_t) +
2757 ck.Advection_strong(dmom_adv_sge,grad_u) +
2758 ck.Hamiltonian_strong(dmom_ham_grad_sge,grad_u) +
2759 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
2760 ck.Reaction_strong(mom_u_source) -
2761 ck.Reaction_strong(dmom_u_acc_u*
u*div_mesh_velocity);
2763 pdeResidual_v =
ck.Mass_strong(mom_v_acc_t) +
2764 ck.Advection_strong(dmom_adv_sge,grad_v) +
2765 ck.Hamiltonian_strong(dmom_ham_grad_sge,grad_v) +
2766 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
2767 ck.Reaction_strong(mom_v_source) -
2768 ck.Reaction_strong(dmom_v_acc_v*
v*div_mesh_velocity);
2772 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
2774 elementDiameter.data()[eN],
2779 dmom_u_ham_grad_p[0],
2782 q_cfl.data()[eN_k]);
2789 dmom_u_ham_grad_p[0],
2792 q_cfl.data()[eN_k]);
2794 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
2795 tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0;
2808 dmom_adv_star[0] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
2809 dmom_adv_star[1] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
2811 mom_u_adv[0] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+0]);
2812 mom_u_adv[1] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+0]);
2814 mom_v_adv[0] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+1]);
2815 mom_v_adv[1] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+1]);
2818 for (
int i=0;i<nDOF_test_element;i++)
2820 int i_nSpace = i*nSpace;
2821 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
2822 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
2824 for (
int i=0;i<nDOF_v_test_element;i++)
2826 int i_nSpace = i*nSpace;
2828 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2829 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2830 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2831 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2834 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2835 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2839 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v);
2840 q_numDiff_u.data()[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2841 (1.0-useMetrics)*hFactor*hFactor*elementDiameter.data()[eN]*elementDiameter.data()[eN]);
2842 q_numDiff_v.data()[eN_k] = q_numDiff_u.data()[eN_k];
2843 q_numDiff_w.data()[eN_k] = q_numDiff_u.data()[eN_k];
2844 numDiffMax = std::fmax(q_numDiff_u.data()[eN_k], numDiffMax);
2848 double level_set_normal[nSpace];
2852 double norm_exact=0.0,norm_cut=0.0;
2853 if (use_ball_as_particle)
2855 for (
int I=0;I<nSpace;I++)
2859 norm_cut += level_set_normal[I]*level_set_normal[I];
2860 norm_exact += ball_n[I]*ball_n[I];
2865 for (
int I=0;I<nSpace;I++)
2869 norm_cut += level_set_normal[I]*level_set_normal[I];
2870 norm_exact += particle_signed_distance_normals.data()[eN_k_3d+I]*particle_signed_distance_normals.data()[eN_k_3d+I];
2873 norm_cut = std::sqrt(norm_cut);
2874 norm_exact = std::sqrt(norm_exact);
2875 assert(std::fabs(1.0-norm_cut) < 1.0e-8);
2876 assert(std::fabs(1.0-norm_exact) < 1.0e-8);
2878 for (
int I=0;I<nSpace;I++)
2879 level_set_normal[I]*=-1.0;
2889 if (use_ball_as_particle)
2890 for (
int I=0;I<nSpace;I++)
2891 level_set_normal[I] = ball_n[I];
2893 for (
int I=0;I<nSpace;I++)
2894 level_set_normal[I] = particle_signed_distance_normals.data()[eN_k_3d+I];
2897 NONCONSERVATIVE_FORM,
2898 eN < nElements_owned,
2902 nQuadraturePoints_global,
2903 &particle_signed_distances.data()[eN_k],
2905 &particle_velocities.data()[eN_k_3d],
2906 particle_centroids.data(),
2907 use_ball_as_particle,
2910 ball_velocity.data(),
2911 ball_angular_velocity.data(),
2912 ball_density.data(),
2914 particle_penalty_constant/h_phi,
2933 q_velocity_sge.data()[eN_k_nSpace+0],
2934 q_velocity_sge.data()[eN_k_nSpace+1],
2935 q_velocity_sge.data()[eN_k_nSpace+1],
2954 dmom_u_ham_grad_u_s,
2955 dmom_u_ham_grad_v_s,
2960 dmom_v_ham_grad_u_s,
2961 dmom_v_ham_grad_v_s,
2966 dmom_w_ham_grad_w_s,
2974 particle_netForces.data(),
2975 particle_netMoments.data(),
2976 particle_surfaceArea.data(),
2977 particle_surfaceArea_projected.data(),
2978 projection_direction.data(),
2979 particle_volume.data());
2995 q_mom_u_acc.data()[eN_k] = mom_u_acc;
2996 q_mom_v_acc.data()[eN_k] = mom_v_acc;
2998 q_mass_adv.data()[eN_k_nSpace+0] =
u;
2999 q_mass_adv.data()[eN_k_nSpace+1] =
v;
3003 if (use_ball_as_particle)
3005 q_mom_u_acc.data()[eN_k] = particle_velocities.data()[eN_k_3d+0];
3006 q_mom_v_acc.data()[eN_k] = particle_velocities.data()[eN_k_3d+1];
3007 q_mass_adv.data()[eN_k_nSpace+0] = particle_velocities.data()[eN_k_3d+0];
3008 q_mass_adv.data()[eN_k_nSpace+1] = particle_velocities.data()[eN_k_3d+1];
3012 q_mom_u_acc.data()[eN_k] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+0];
3013 q_mom_v_acc.data()[eN_k] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+1];
3014 q_mass_adv.data()[eN_k_nSpace+0] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+0];
3015 q_mass_adv.data()[eN_k_nSpace+1] = particle_velocities.data()[particle_index*nQuadraturePoints_global + eN_k_3d+1];
3021 double mesh_vel[nSpace];
3027 if (fluid_phase == 0)
3028 H_f =
gf.
ImH(0.,0.);
3042 if ((eN < nElements_owned) && isActiveElement[eN])
3044 domain_volume += H_s*dV*H_f;
3045 p_L1 += fabs(p_e)*H_s*dV*H_f;
3046 u_L1 += fabs(u_e)*H_s*dV*H_f;
3047 v_L1 += fabs(v_e)*H_s*dV*H_f;
3048 velocity_L1 += fabs(velocity_e)*H_s*dV*H_f;
3050 p_L2 += p_e*p_e*H_s*dV*H_f;
3051 u_L2 += u_e*u_e*H_s*dV*H_f;
3052 v_L2 += v_e*v_e*H_s*dV*H_f;
3053 velocity_L2 += velocity_e*velocity_e*H_s*dV*H_f;
3054 p_dv += p*H_s*H_f*dV;
3055 pa_dv += q_u_0.data()[eN_k]*H_s*H_f*dV;
3056 total_volume+=H_s*H_f*dV;
3057 total_surface_area+=D_s*H_f*dV;
3058 if (phi_solid.data()[eN_k] >= 0.0)
3060 p_LI = fmax(p_LI, fabs(p_e));
3061 u_LI = fmax(u_LI, fabs(u_e));
3062 v_LI = fmax(v_LI, fabs(v_e));
3063 velocity_LI = fmax(velocity_LI, fabs(velocity_e));
3066 for(
int i=0;i<nDOF_test_element;i++)
3068 int i_nSpace=i*nSpace;
3069 elementResidual_mesh[i] += H_s*H_f*(
ck.Reaction_weak(1.0,p_test_dV[i]) -
3070 ck.Reaction_weak(1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) -
3071 ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace]));
3072 elementResidual_p[i] += H_s*H_f*(
ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace])
3073 +
ck.Hamiltonian_weak(mass_ham, p_test_dV[i])
3074 +
DM*MOVING_DOMAIN*(
ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]) -
3075 ck.Reaction_weak(alphaBDF*1.0,p_test_dV[i]*q_dV_last.data()[eN_k]/dV) -
3076 ck.Advection_weak(mesh_vel,&p_grad_test_dV[i_nSpace])) +
3077 ck.Reaction_weak(mass_source,p_test_dV[i]));
3078 if (nDOF_test_element == nDOF_v_test_element)
3080 elementResidual_p[i] +=
3081 H_s*H_f*(PRESSURE_PROJECTION_STABILIZATION *
ck.pressureProjection_weak(mom_uu_diff_ten[1], p, p_element_avg, p_test_ref.data()[k*nDOF_test_element+i], dV) +
3082 (1 - PRESSURE_PROJECTION_STABILIZATION) *
ck.SubgridError(subgridError_u,Lstar_u_p[i]) +
3083 (1 - PRESSURE_PROJECTION_STABILIZATION) *
ck.SubgridError(subgridError_v,Lstar_v_p[i]));
3085 if (PRESSURE_PROJECTION_STABILIZATION==1. && mom_uu_diff_ten[1]==0.)
3087 printf(
"Warning the Bochev-Dohrnmann-Gunzburger stabilization cannot be applied to inviscid fluids.");
3091 if (
gf_s.
D(0.,0.) == 0.0)
3092 assert(mass_source_s == 0.0);
3093 elementResidual_p[i] += H_f*(
ck.Reaction_weak(mass_source_s,p_test_dV[i]));
3096 for(
int i=0;i<nDOF_v_test_element;i++)
3098 int i_nSpace=i*nSpace;
3099 elementResidual_u[i] += H_s*H_f*(
ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
3100 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
3101 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]) +
3102 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]) +
3103 ck.Reaction_weak(mom_u_source+NONCONSERVATIVE_FORM*dmom_u_acc_u*
u*div_mesh_velocity,vel_test_dV[i]) +
3104 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
3105 MOMENTUM_SGE*VELOCITY_SGE*
ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
3106 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]));
3107 elementResidual_v[i] += H_s*H_f*(
ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
3108 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
3109 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]) +
3110 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]) +
3111 ck.Reaction_weak(mom_v_source+NONCONSERVATIVE_FORM*dmom_v_acc_v*
v*div_mesh_velocity,vel_test_dV[i]) +
3112 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
3113 MOMENTUM_SGE*VELOCITY_SGE*
ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
3114 ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]));
3115 elementResidual_u[i] += H_s*H_f*MOMENTUM_SGE*PRESSURE_SGE*
ck.SubgridError(subgridError_p,Lstar_p_u[i]);
3116 elementResidual_v[i] += H_s*H_f*MOMENTUM_SGE*PRESSURE_SGE*
ck.SubgridError(subgridError_p,Lstar_p_v[i]);
3119 elementResidual_u[i] += H_f*(
ck.Advection_weak(mom_u_adv_s,&vel_grad_test_dV[i_nSpace]) +
3120 ck.Reaction_weak(mom_u_source_s,vel_test_dV[i]) +
3121 ck.Hamiltonian_weak(mom_u_ham_s,vel_test_dV[i]));
3122 elementResidual_v[i] += H_f*(
ck.Advection_weak(mom_v_adv_s,&vel_grad_test_dV[i_nSpace]) +
3123 ck.Reaction_weak(mom_v_source_s,vel_test_dV[i]) +
3124 ck.Hamiltonian_weak(mom_v_ham_s,vel_test_dV[i]));
3128 numerical_viscosity.data()[eN_k] = q_numDiff_u_last.data()[eN_k] + MOMENTUM_SGE*VELOCITY_SGE*tau_v*(dmom_adv_star[0]*dmom_adv_star[0]+
3129 dmom_adv_star[1]*dmom_adv_star[1]);
3130 if (!isActiveElement[eN])
3132 assert(std::fabs(
gf_s.
H(particle_eps,phi_solid.data()[eN_k])) == 0.0);
3133 assert(std::fabs(
gf_s.
D(particle_eps,phi_solid.data()[eN_k])) == 0.0);
3138 for(
int k=0;k<nQuadraturePoints_element;k++)
3141 int eN_k = eN*nQuadraturePoints_element+k;
3142 q_numDiff_u.data()[eN_k] = numDiffMax;
3143 q_numDiff_v.data()[eN_k] = numDiffMax;
3144 q_numDiff_w.data()[eN_k] = numDiffMax;
3150 for(
int i=0;i<nDOF_test_element;i++)
3152 int eN_i=eN*nDOF_test_element+i;
3153 elementResidual_p_save.data()[eN_i] += elementResidual_p[i];
3154 mesh_volume_conservation_element_weak += elementResidual_mesh[i];
3155 if (!isActiveElement[eN])
3157 assert(elementResidual_p[i]==0.0);
3159 globalResidual.data()[offset_p+stride_p*rp_l2g.data()[eN_i]]+=elementResidual_p[i];
3162 isActiveR.data()[offset_p+stride_p*rp_l2g.data()[eN_i]] = 1.0;
3163 isActiveDOF_p.data()[p_l2g.data()[eN_i]] = 1.0;
3166 for(
int i=0;i<nDOF_v_test_element;i++)
3168 int eN_i=eN*nDOF_v_test_element+i;
3169 if (!isActiveElement[eN])
3171 assert(elementResidual_u[i]==0.0);
3172 assert(elementResidual_v[i]==0.0);
3174 globalResidual.data()[offset_u+stride_u*rvel_l2g.data()[eN_i]]+=elementResidual_u[i];
3175 globalResidual.data()[offset_v+stride_v*rvel_l2g.data()[eN_i]]+=elementResidual_v[i];
3178 isActiveR.data()[offset_u+stride_u*rvel_l2g.data()[eN_i]] = 1.0;
3179 isActiveR.data()[offset_v+stride_v*rvel_l2g.data()[eN_i]] = 1.0;
3180 isActiveDOF_vel.data()[vel_l2g.data()[eN_i]] = 1.0;
3182 double x = mesh_dof.data()[3*mesh_l2g.data()[eN_i]+0],
3183 y = mesh_dof.data()[3*mesh_l2g.data()[eN_i]+1],
3184 z = mesh_dof.data()[3*mesh_l2g.data()[eN_i]+2];
3187 ball_velocity.data(),ball_angular_velocity.data(),
3188 particle_index,x,y,
z,
3189 ball_u.data()[vel_l2g.data()[eN_i]],ball_v.data()[vel_l2g.data()[eN_i]]);
3191 mesh_volume_conservation += mesh_volume_conservation_element;
3192 mesh_volume_conservation_weak += mesh_volume_conservation_element_weak;
3193 mesh_volume_conservation_err_max=fmax(mesh_volume_conservation_err_max,fabs(mesh_volume_conservation_element));
3194 mesh_volume_conservation_err_max_weak=fmax(mesh_volume_conservation_err_max_weak,fabs(mesh_volume_conservation_element_weak));
3199 if(isActiveElement[elementBoundaryElementsArray[(*it)*2+0]] && isActiveElement[elementBoundaryElementsArray[(*it)*2+1]])
3201 std::map<int,double> DWp_Dn_jump, DW_Dn_jump;
3202 double gamma_cutfem=ghost_penalty_constant,gamma_cutfem_p=ghost_penalty_constant,h_cutfem=elementBoundaryDiameter.data()[*it];
3203 int eN_nDOF_v_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_v_trial_element;
3207 for (
int i_offset=1;i_offset<nDOF_v_trial_element;i_offset++)
3210 double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
3211 v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
3212 norm_v=fmax(norm_v,sqrt(
u*
u+
v*
v));
3214 double gamma_v_dim =
rho_0*(
nu_0 + norm_v*h_cutfem + alphaBDF*h_cutfem*h_cutfem);
3215 gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
3216 if (NONCONSERVATIVE_FORM)
3217 gamma_cutfem*=gamma_v_dim;
3219 gamma_cutfem*=(gamma_v_dim/
rho_0);
3220 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3222 double Dp_Dn_jump=0.0, Du_Dn_jump=0.0, Dv_Dn_jump=0.0,dS;
3223 for (
int eN_side=0;eN_side < 2; eN_side++)
3226 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
3227 for (
int i=0;i<nDOF_test_element;i++)
3229 DWp_Dn_jump[rp_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
3231 for (
int i=0;i<nDOF_v_test_element;i++)
3233 DW_Dn_jump[rvel_l2g.data()[eN*nDOF_v_test_element+i]] = 0.0;
3236 for (
int eN_side=0;eN_side < 2; eN_side++)
3239 eN = elementBoundaryElementsArray[ebN*2+eN_side],
3240 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+eN_side],
3241 eN_nDOF_trial_element = eN*nDOF_trial_element,
3242 eN_nDOF_v_trial_element = eN*nDOF_v_trial_element,
3243 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3244 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3251 jac_int[nSpace*nSpace],
3253 jacInv_int[nSpace*nSpace],
3254 boundaryJac[nSpace*(nSpace-1)],
3255 metricTensor[(nSpace-1)*(nSpace-1)],
3256 metricTensorDetSqrt,
3257 p_test_dS[nDOF_test_element],vel_test_dS[nDOF_v_test_element],
3258 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_v_trial_element*nSpace],
3259 p_grad_test_dS[nDOF_trial_element*nSpace],vel_grad_test_dS[nDOF_v_trial_element*nSpace],
3260 normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
3261 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3262 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;
3264 ck.calculateMapping_elementBoundary(eN,
3270 mesh_trial_trace_ref.data(),
3271 mesh_grad_trial_trace_ref.data(),
3272 boundaryJac_ref.data(),
3278 metricTensorDetSqrt,
3283 ck.calculateMappingVelocity_elementBoundary(eN,
3287 mesh_velocity_dof.data(),
3289 mesh_trial_trace_ref.data(),
3290 xt_int,yt_int,zt_int,
3295 dS = metricTensorDetSqrt*dS_ref.data()[kb];
3298 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,p_grad_trial_trace);
3299 ck_v.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_v_trial_element],jacInv_int,vel_grad_trial_trace);
3301 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_int);
3302 ck_v.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],u_int);
3303 ck_v.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],v_int);
3304 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_int);
3305 ck_v.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_int);
3306 ck_v.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_int);
3307 for (
int I=0;I<nSpace;I++)
3309 Dp_Dn_jump += grad_p_int[I]*normal[I];
3310 Du_Dn_jump += grad_u_int[I]*normal[I];
3311 Dv_Dn_jump += grad_v_int[I]*normal[I];
3313 for (
int i=0;i<nDOF_test_element;i++)
3315 for (
int I=0;I<nSpace;I++)
3316 DWp_Dn_jump[rp_l2g[eN_nDOF_trial_element+i]] += p_grad_trial_trace[i*nSpace+I]*normal[I];
3318 for (
int i=0;i<nDOF_v_test_element;i++)
3320 for (
int I=0;I<nSpace;I++)
3321 DW_Dn_jump[rvel_l2g[eN_nDOF_v_trial_element+i]] += vel_grad_trial_trace[i*nSpace+I]*normal[I];
3324 for (std::map<int,double>::iterator W_it=DWp_Dn_jump.begin(); W_it!=DWp_Dn_jump.end(); ++W_it)
3326 int i_global = W_it->first;
3327 double DWp_Dn_jump_i = W_it->second;
3328 globalResidual.data()[offset_p+stride_p*i_global]+=gamma_cutfem_p*h_cutfem*Dp_Dn_jump*DWp_Dn_jump_i*dS;
3330 for (std::map<int,double>::iterator W_it=DW_Dn_jump.begin(); W_it!=DW_Dn_jump.end(); ++W_it)
3332 int i_global = W_it->first;
3333 double DW_Dn_jump_i = W_it->second;
3334 globalResidual.data()[offset_u+stride_u*i_global]+=gamma_cutfem*h_cutfem*Du_Dn_jump*DW_Dn_jump_i*dS;
3335 globalResidual.data()[offset_v+stride_v*i_global]+=gamma_cutfem*h_cutfem*Dv_Dn_jump*DW_Dn_jump_i*dS;
3351 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3353 int ebN = exteriorElementBoundariesArray.data()[ebNE],
3354 eN = elementBoundaryElementsArray.data()[ebN*2+0],
3355 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
3356 eN_nDOF_trial_element = eN*nDOF_trial_element,
3357 eN_nDOF_v_trial_element = eN*nDOF_v_trial_element;
3358 if (boundaryFlags[ebN] < 1)
3360 double elementResidual_mesh[nDOF_test_element],
3361 elementResidual_p[nDOF_test_element],
3362 elementResidual_u[nDOF_v_test_element],
3363 elementResidual_v[nDOF_v_test_element],
3365 const double* elementResidual_w(NULL);
3366 for (
int i=0;i<nDOF_test_element;i++)
3368 elementResidual_mesh[i]=0.0;
3369 elementResidual_p[i]=0.0;
3371 for (
int i=0;i<nDOF_v_test_element;i++)
3373 elementResidual_u[i]=0.0;
3374 elementResidual_v[i]=0.0;
3376 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
3377 for (
int j=0;j<nDOF_mesh_trial_element;j++)
3379 int eN_j = eN*nDOF_mesh_trial_element+j;
3380 element_phi[j] = phi_nodes.data()[p_l2g.data()[eN_j]];
3381 element_phi_s[j] = phi_solid_nodes[p_l2g.data()[eN_j]];
3383 double element_nodes[nDOF_mesh_trial_element*3];
3384 for (
int i=0;i<nDOF_mesh_trial_element;i++)
3386 int eN_i=eN*nDOF_mesh_trial_element+i;
3387 for(
int I=0;I<3;I++)
3388 element_nodes[i*3 + I] = mesh_dof[mesh_l2g.data()[eN_i]*3 + I];
3390 double mesh_dof_ref[nDOF_mesh_trial_element*3]={0.,0.,0.,1.,0.,0.,0.,1.,0.};
3391 double xb_ref_calc[nQuadraturePoints_elementBoundary*3];
3392 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3394 double x=0.0,y=0.0,
z=0.0;
3395 for (
int j=0;j<nDOF_mesh_trial_element;j++)
3397 int ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb;
3398 int ebN_local_kb_j = ebN_local_kb*nDOF_mesh_trial_element+j;
3399 x += mesh_dof_ref[j*3+0]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
3400 y += mesh_dof_ref[j*3+1]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
3401 z += mesh_dof_ref[j*3+2]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
3403 xb_ref_calc[3*kb+0] = x;
3404 xb_ref_calc[3*kb+1] = y;
3405 xb_ref_calc[3*kb+2] =
z;
3407 int icase_s =
gf_s.
calculate(element_phi_s, element_nodes, xb_ref_calc,
true);
3409 int icase =
gf.
calculate(element_phi, element_nodes, xb_ref.data(), -
rho_1*g.data()[1], -
rho_0*g.data()[1],
true,
true);
3411 int icase =
gf.
calculate(element_phi, element_nodes, xb_ref.data(), 1.0,1.0,
true,
false);
3414 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3416 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3417 ebNE_kb_nSpace = ebNE_kb*nSpace,
3418 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3419 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3420 double phi_s_ext=0.0,
3429 p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0,
3432 dmom_u_acc_u_ext=0.0,
3434 dmom_v_acc_v_ext=0.0,
3436 dmom_w_acc_w_ext=0.0,
3438 dmass_adv_u_ext[nSpace]=
ZEROVEC,
3439 dmass_adv_v_ext[nSpace]=
ZEROVEC,
3440 dmass_adv_w_ext[nSpace]=
ZEROVEC,
3441 mom_u_adv_ext[nSpace]=
ZEROVEC,
3442 dmom_u_adv_u_ext[nSpace]=
ZEROVEC,
3443 dmom_u_adv_v_ext[nSpace]=
ZEROVEC,
3444 dmom_u_adv_w_ext[nSpace]=
ZEROVEC,
3445 mom_v_adv_ext[nSpace]=
ZEROVEC,
3446 dmom_v_adv_u_ext[nSpace]=
ZEROVEC,
3447 dmom_v_adv_v_ext[nSpace]=
ZEROVEC,
3448 dmom_v_adv_w_ext[nSpace]=
ZEROVEC,
3449 mom_w_adv_ext[nSpace]=
ZEROVEC,
3450 dmom_w_adv_u_ext[nSpace]=
ZEROVEC,
3451 dmom_w_adv_v_ext[nSpace]=
ZEROVEC,
3452 dmom_w_adv_w_ext[nSpace]=
ZEROVEC,
3453 mom_uu_diff_ten_ext[nSpace]=
ZEROVEC,
3454 mom_vv_diff_ten_ext[nSpace]=
ZEROVEC,
3455 mom_ww_diff_ten_ext[nSpace]=
ZEROVEC,
3456 mom_uv_diff_ten_ext[1],
3457 mom_uw_diff_ten_ext[1],
3458 mom_vu_diff_ten_ext[1],
3459 mom_vw_diff_ten_ext[1],
3460 mom_wu_diff_ten_ext[1],
3461 mom_wv_diff_ten_ext[1],
3462 mom_u_source_ext=0.0,
3463 mom_v_source_ext=0.0,
3464 mom_w_source_ext=0.0,
3466 dmom_u_ham_grad_p_ext[nSpace]=
ZEROVEC,
3467 dmom_u_ham_grad_u_ext[nSpace]=
ZEROVEC,
3468 dmom_u_ham_u_ext=0.0,
3469 dmom_u_ham_v_ext=0.0,
3470 dmom_u_ham_w_ext=0.0,
3472 dmom_v_ham_grad_p_ext[nSpace]=
ZEROVEC,
3473 dmom_v_ham_grad_v_ext[nSpace]=
ZEROVEC,
3474 dmom_v_ham_u_ext=0.0,
3475 dmom_v_ham_v_ext=0.0,
3476 dmom_v_ham_w_ext=0.0,
3478 dmom_w_ham_grad_p_ext[nSpace]=
ZEROVEC,
3479 dmom_w_ham_grad_w_ext[nSpace]=
ZEROVEC,
3480 dmom_w_ham_u_ext=0.0,
3481 dmom_w_ham_v_ext=0.0,
3482 dmom_w_ham_w_ext=0.0,
3483 dmom_u_adv_p_ext[nSpace]=
ZEROVEC,
3484 dmom_v_adv_p_ext[nSpace]=
ZEROVEC,
3485 dmom_w_adv_p_ext[nSpace]=
ZEROVEC,
3487 flux_mom_u_adv_ext=0.0,
3488 flux_mom_v_adv_ext=0.0,
3489 flux_mom_w_adv_ext=0.0,
3490 flux_mom_uu_diff_ext=0.0,
3491 flux_mom_uv_diff_ext=0.0,
3492 flux_mom_uw_diff_ext=0.0,
3493 flux_mom_vu_diff_ext=0.0,
3494 flux_mom_vv_diff_ext=0.0,
3495 flux_mom_vw_diff_ext=0.0,
3496 flux_mom_wu_diff_ext=0.0,
3497 flux_mom_wv_diff_ext=0.0,
3498 flux_mom_ww_diff_ext=0.0,
3503 bc_mom_u_acc_ext=0.0,
3504 bc_dmom_u_acc_u_ext=0.0,
3505 bc_mom_v_acc_ext=0.0,
3506 bc_dmom_v_acc_v_ext=0.0,
3507 bc_mom_w_acc_ext=0.0,
3508 bc_dmom_w_acc_w_ext=0.0,
3509 bc_mass_adv_ext[nSpace]=
ZEROVEC,
3510 bc_dmass_adv_u_ext[nSpace]=
ZEROVEC,
3511 bc_dmass_adv_v_ext[nSpace]=
ZEROVEC,
3512 bc_dmass_adv_w_ext[nSpace]=
ZEROVEC,
3513 bc_mom_u_adv_ext[nSpace]=
ZEROVEC,
3514 bc_dmom_u_adv_u_ext[nSpace]=
ZEROVEC,
3515 bc_dmom_u_adv_v_ext[nSpace]=
ZEROVEC,
3516 bc_dmom_u_adv_w_ext[nSpace]=
ZEROVEC,
3517 bc_mom_v_adv_ext[nSpace]=
ZEROVEC,
3518 bc_dmom_v_adv_u_ext[nSpace]=
ZEROVEC,
3519 bc_dmom_v_adv_v_ext[nSpace]=
ZEROVEC,
3520 bc_dmom_v_adv_w_ext[nSpace]=
ZEROVEC,
3521 bc_mom_w_adv_ext[nSpace]=
ZEROVEC,
3522 bc_dmom_w_adv_u_ext[nSpace]=
ZEROVEC,
3523 bc_dmom_w_adv_v_ext[nSpace]=
ZEROVEC,
3524 bc_dmom_w_adv_w_ext[nSpace]=
ZEROVEC,
3525 bc_mom_uu_diff_ten_ext[nSpace]=
ZEROVEC,
3526 bc_mom_vv_diff_ten_ext[nSpace]=
ZEROVEC,
3527 bc_mom_ww_diff_ten_ext[nSpace]=
ZEROVEC,
3528 bc_mom_uv_diff_ten_ext[1],
3529 bc_mom_uw_diff_ten_ext[1],
3530 bc_mom_vu_diff_ten_ext[1],
3531 bc_mom_vw_diff_ten_ext[1],
3532 bc_mom_wu_diff_ten_ext[1],
3533 bc_mom_wv_diff_ten_ext[1],
3534 bc_mom_u_source_ext=0.0,
3535 bc_mom_v_source_ext=0.0,
3536 bc_mom_w_source_ext=0.0,
3537 bc_mom_u_ham_ext=0.0,
3538 bc_dmom_u_ham_grad_p_ext[nSpace]=
ZEROVEC,
3539 bc_dmom_u_ham_grad_u_ext[nSpace]=
ZEROVEC,
3540 bc_dmom_u_ham_u_ext=0.0,
3541 bc_dmom_u_ham_v_ext=0.0,
3542 bc_dmom_u_ham_w_ext=0.0,
3543 bc_mom_v_ham_ext=0.0,
3544 bc_dmom_v_ham_grad_p_ext[nSpace]=
ZEROVEC,
3545 bc_dmom_v_ham_grad_v_ext[nSpace]=
ZEROVEC,
3546 bc_dmom_v_ham_u_ext=0.0,
3547 bc_dmom_v_ham_v_ext=0.0,
3548 bc_dmom_v_ham_w_ext=0.0,
3549 bc_mom_w_ham_ext=0.0,
3550 bc_dmom_w_ham_grad_p_ext[nSpace]=
ZEROVEC,
3551 bc_dmom_w_ham_grad_w_ext[nSpace]=
ZEROVEC,
3552 bc_dmom_w_ham_u_ext=0.0,
3553 bc_dmom_w_ham_v_ext=0.0,
3554 bc_dmom_w_ham_w_ext=0.0,
3555 jac_ext[nSpace*nSpace],
3557 jacInv_ext[nSpace*nSpace],
3558 boundaryJac[nSpace*(nSpace-1)],
3559 metricTensor[(nSpace-1)*(nSpace-1)],
3560 metricTensorDetSqrt,
3561 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_v_test_element],
3562 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_v_trial_element*nSpace],
3563 vel_grad_test_dS[nDOF_v_trial_element*nSpace],
3564 normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3568 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3569 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;
3573 ck.calculateMapping_elementBoundary(eN,
3579 mesh_trial_trace_ref.data(),
3580 mesh_grad_trial_trace_ref.data(),
3581 boundaryJac_ref.data(),
3587 metricTensorDetSqrt,
3591 ck.calculateMappingVelocity_elementBoundary(eN,
3595 mesh_velocity_dof.data(),
3597 mesh_trial_trace_ref.data(),
3598 xt_ext,yt_ext,zt_ext,
3610 dS = metricTensorDetSqrt*dS_ref.data()[kb];
3613 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3614 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
3616 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3617 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3621 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
3622 ck_v.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_v_trial_element],jacInv_ext,vel_grad_trial_trace);
3624 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext);
3625 ck_v.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],u_ext);
3626 ck_v.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],v_ext);
3627 ck.valFromDOF(p_old_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_old);
3628 ck_v.valFromDOF(u_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],u_old);
3629 ck_v.valFromDOF(v_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],v_old);
3630 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
3631 ck_v.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_ext);
3632 ck_v.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_ext);
3633 ck.gradFromDOF(p_old_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_old);
3634 ck_v.gradFromDOF(u_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_old);
3635 ck_v.gradFromDOF(v_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_old);
3636 ck.valFromDOF(phi_solid_nodes.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],phi_s_ext);
3638 for (
int j=0;j<nDOF_test_element;j++)
3640 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
3642 for (
int j=0;j<nDOF_v_test_element;j++)
3644 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_v_test_element+j]*dS;
3645 for (
int I=0;I<nSpace;I++)
3646 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3648 bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
3650 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*(ebqe_bc_u_ext.data()[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
3651 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*(ebqe_bc_v_ext.data()[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
3653 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
3657 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);
3658 if (use_ball_as_particle == 1 && nParticles > 0)
3660 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x_ext,y_ext,z_ext,ebqe_phi_s.data()[ebNE_kb]);
3663 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3666 double H = (1.0-useVF)*
gf.
H(eps_rho,ebqe_phi_ext[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,ebqe_vf_ext[ebNE_kb]));
3667 double ImH = (1.0-useVF)*
gf.
ImH(eps_rho,ebqe_phi_ext[ebNE_kb]) + useVF*(1.0-fmin(1.0,fmax(0.0,ebqe_vf_ext[ebNE_kb])));
3675 elementDiameter.data()[eN],
3676 smagorinskyConstant,
3677 turbulenceClosureModel,
3680 ebqe_vf_ext.data()[ebNE_kb],
3681 ebqe_phi_ext.data()[ebNE_kb],
3682 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
3683 ebqe_kappa_phi_ext.data()[ebNE_kb],
3687 ebqe_phi_s.data()[ebNE_kb],
3705 ebqe_eddy_viscosity.data()[ebNE_kb],
3706 ebqe_eddy_viscosity_last.data()[ebNE_kb],
3729 mom_uu_diff_ten_ext,
3730 mom_vv_diff_ten_ext,
3731 mom_ww_diff_ten_ext,
3732 mom_uv_diff_ten_ext,
3733 mom_uw_diff_ten_ext,
3734 mom_vu_diff_ten_ext,
3735 mom_vw_diff_ten_ext,
3736 mom_wu_diff_ten_ext,
3737 mom_wv_diff_ten_ext,
3742 dmom_u_ham_grad_p_ext,
3743 dmom_u_ham_grad_u_ext,
3748 dmom_v_ham_grad_p_ext,
3749 dmom_v_ham_grad_v_ext,
3754 dmom_w_ham_grad_p_ext,
3755 dmom_w_ham_grad_w_ext,
3763 H = (1.0-useVF)*
gf.
H(eps_rho,bc_ebqe_phi_ext[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,bc_ebqe_vf_ext[ebNE_kb]));
3764 ImH = (1.0-useVF)*
gf.
ImH(eps_rho,bc_ebqe_phi_ext[ebNE_kb]) + useVF*(1.0-fmin(1.0,fmax(0.0,bc_ebqe_vf_ext[ebNE_kb])));
3772 elementDiameter.data()[eN],
3773 smagorinskyConstant,
3774 turbulenceClosureModel,
3777 bc_ebqe_vf_ext.data()[ebNE_kb],
3778 bc_ebqe_phi_ext.data()[ebNE_kb],
3779 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
3780 ebqe_kappa_phi_ext.data()[ebNE_kb],
3784 ebqe_phi_s.data()[ebNE_kb],
3802 bc_eddy_viscosity_ext,
3803 ebqe_eddy_viscosity_last.data()[ebNE_kb],
3805 bc_dmom_u_acc_u_ext,
3807 bc_dmom_v_acc_v_ext,
3809 bc_dmom_w_acc_w_ext,
3815 bc_dmom_u_adv_u_ext,
3816 bc_dmom_u_adv_v_ext,
3817 bc_dmom_u_adv_w_ext,
3819 bc_dmom_v_adv_u_ext,
3820 bc_dmom_v_adv_v_ext,
3821 bc_dmom_v_adv_w_ext,
3823 bc_dmom_w_adv_u_ext,
3824 bc_dmom_w_adv_v_ext,
3825 bc_dmom_w_adv_w_ext,
3826 bc_mom_uu_diff_ten_ext,
3827 bc_mom_vv_diff_ten_ext,
3828 bc_mom_ww_diff_ten_ext,
3829 bc_mom_uv_diff_ten_ext,
3830 bc_mom_uw_diff_ten_ext,
3831 bc_mom_vu_diff_ten_ext,
3832 bc_mom_vw_diff_ten_ext,
3833 bc_mom_wu_diff_ten_ext,
3834 bc_mom_wv_diff_ten_ext,
3835 bc_mom_u_source_ext,
3836 bc_mom_v_source_ext,
3837 bc_mom_w_source_ext,
3839 bc_dmom_u_ham_grad_p_ext,
3840 bc_dmom_u_ham_grad_u_ext,
3841 bc_dmom_u_ham_u_ext,
3842 bc_dmom_u_ham_v_ext,
3843 bc_dmom_u_ham_w_ext,
3845 bc_dmom_v_ham_grad_p_ext,
3846 bc_dmom_v_ham_grad_v_ext,
3847 bc_dmom_v_ham_u_ext,
3848 bc_dmom_v_ham_v_ext,
3849 bc_dmom_v_ham_w_ext,
3851 bc_dmom_w_ham_grad_p_ext,
3852 bc_dmom_w_ham_grad_w_ext,
3853 bc_dmom_w_ham_u_ext,
3854 bc_dmom_w_ham_v_ext,
3855 bc_dmom_w_ham_w_ext,
3861 if (turbulenceClosureModel >= 3)
3863 const double turb_var_grad_0_dummy[nSpace] =
ZEROVEC;
3864 const double c_mu = 0.09;
3866 turbulenceClosureModel,
3874 ebqe_vf_ext.data()[ebNE_kb],
3875 ebqe_phi_ext.data()[ebNE_kb],
3878 ebqe_turb_var_0.data()[ebNE_kb],
3879 ebqe_turb_var_1.data()[ebNE_kb],
3880 turb_var_grad_0_dummy,
3881 ebqe_eddy_viscosity.data()[ebNE_kb],
3882 mom_uu_diff_ten_ext,
3883 mom_vv_diff_ten_ext,
3884 mom_ww_diff_ten_ext,
3885 mom_uv_diff_ten_ext,
3886 mom_uw_diff_ten_ext,
3887 mom_vu_diff_ten_ext,
3888 mom_vw_diff_ten_ext,
3889 mom_wu_diff_ten_ext,
3890 mom_wv_diff_ten_ext,
3896 turbulenceClosureModel,
3904 bc_ebqe_vf_ext.data()[ebNE_kb],
3905 bc_ebqe_phi_ext.data()[ebNE_kb],
3908 ebqe_turb_var_0.data()[ebNE_kb],
3909 ebqe_turb_var_1.data()[ebNE_kb],
3910 turb_var_grad_0_dummy,
3911 bc_eddy_viscosity_ext,
3912 bc_mom_uu_diff_ten_ext,
3913 bc_mom_vv_diff_ten_ext,
3914 bc_mom_ww_diff_ten_ext,
3915 bc_mom_uv_diff_ten_ext,
3916 bc_mom_uw_diff_ten_ext,
3917 bc_mom_vu_diff_ten_ext,
3918 bc_mom_vw_diff_ten_ext,
3919 bc_mom_wu_diff_ten_ext,
3920 bc_mom_wv_diff_ten_ext,
3921 bc_mom_u_source_ext,
3922 bc_mom_v_source_ext,
3923 bc_mom_w_source_ext);
3930 if (NONCONSERVATIVE_FORM > 0.0)
3932 mom_u_ham_ext -= MOVING_DOMAIN*dmom_u_acc_u_ext*(grad_u_ext[0]*xt_ext + grad_u_ext[1]*yt_ext);
3933 dmom_u_ham_grad_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
3934 dmom_u_ham_grad_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
3938 mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext;
3939 mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext;
3940 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
3941 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
3945 if (NONCONSERVATIVE_FORM > 0.0)
3947 mom_v_ham_ext -= MOVING_DOMAIN*dmom_v_acc_v_ext*(grad_v_ext[0]*xt_ext + grad_v_ext[1]*yt_ext);
3948 dmom_v_ham_grad_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
3949 dmom_v_ham_grad_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
3953 mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext;
3954 mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext;
3955 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
3956 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
3960 if (NONCONSERVATIVE_FORM < 1.0)
3962 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext;
3963 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext;
3965 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext;
3966 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext;
3971 ck.calculateGScale(G,normal,h_penalty);
3972 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
3974 isDOFBoundary_p.data()[ebNE_kb],
3975 isDOFBoundary_u.data()[ebNE_kb],
3976 isDOFBoundary_v.data()[ebNE_kb],
3977 isDOFBoundary_w.data()[ebNE_kb],
3978 isAdvectiveFluxBoundary_p.data()[ebNE_kb],
3979 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
3980 isAdvectiveFluxBoundary_v.data()[ebNE_kb],
3981 isAdvectiveFluxBoundary_w.data()[ebNE_kb],
3982 dmom_u_ham_grad_p_ext[0],
3983 bc_dmom_u_ham_grad_p_ext[0],
3992 ebqe_bc_flux_mass_ext.data()[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),
3993 ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
3994 ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
3995 ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
4007 dmom_u_ham_grad_u_ext,
4023 &ebqe_velocity.data()[ebNE_kb_nSpace]);
4024 for (
int I=0;I<nSpace;I++)
4025 ebqe_velocity.data()[ebNE_kb_nSpace+I]/=porosity_ext;
4027 ebqe_phi_ext.data()[ebNE_kb],
4028 sdInfo_u_u_rowptr.data(),
4029 sdInfo_u_u_colind.data(),
4030 isDOFBoundary_u.data()[ebNE_kb],
4031 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4033 bc_mom_uu_diff_ten_ext,
4035 ebqe_bc_flux_u_diff_ext.data()[ebNE_kb],
4036 mom_uu_diff_ten_ext,
4040 flux_mom_uu_diff_ext);
4042 ebqe_phi_ext.data()[ebNE_kb],
4043 sdInfo_u_v_rowptr.data(),
4044 sdInfo_u_v_colind.data(),
4045 isDOFBoundary_v.data()[ebNE_kb],
4046 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4048 bc_mom_uv_diff_ten_ext,
4051 mom_uv_diff_ten_ext,
4055 flux_mom_uv_diff_ext);
4057 ebqe_phi_ext.data()[ebNE_kb],
4058 sdInfo_v_u_rowptr.data(),
4059 sdInfo_v_u_colind.data(),
4060 isDOFBoundary_u.data()[ebNE_kb],
4061 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4063 bc_mom_vu_diff_ten_ext,
4066 mom_vu_diff_ten_ext,
4070 flux_mom_vu_diff_ext);
4072 ebqe_phi_ext.data()[ebNE_kb],
4073 sdInfo_v_v_rowptr.data(),
4074 sdInfo_v_v_colind.data(),
4075 isDOFBoundary_v.data()[ebNE_kb],
4076 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4078 bc_mom_vv_diff_ten_ext,
4080 ebqe_bc_flux_v_diff_ext.data()[ebNE_kb],
4081 mom_vv_diff_ten_ext,
4085 flux_mom_vv_diff_ext);
4086 flux.data()[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
4094 if (ebN < nElementBoundaries_owned)
4096 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];
4097 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];
4099 force_p_x = p_ext*normal[0];
4100 force_p_y = p_ext*normal[1];
4102 force_x = force_p_x + force_v_x;
4103 force_y = force_p_y + force_v_y;
4105 r_x = x_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+0];
4106 r_y = y_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+1];
4108 wettedAreas.data()[boundaryFlags.data()[ebN]] += dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
4110 netForces_p.data()[3*boundaryFlags.data()[ebN]+0] += force_p_x*dS;
4111 netForces_p.data()[3*boundaryFlags.data()[ebN]+1] += force_p_y*dS;
4113 netForces_v.data()[3*boundaryFlags.data()[ebN]+0] += force_v_x*dS;
4114 netForces_v.data()[3*boundaryFlags.data()[ebN]+1] += force_v_y*dS;
4116 netMoments.data()[3*boundaryFlags.data()[ebN]+2] += (r_x*force_y - r_y*force_x)*dS;
4121 const double H_s =
gf_s.
H(particle_eps, ebqe_phi_s.data()[ebNE_kb]);
4122 if (isActiveElement[eN])
4124 total_flux += flux_mass_ext*dS;
4125 for (
int i=0;i<nDOF_test_element;i++)
4127 elementResidual_mesh[i] -= H_s*
ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),p_test_dS[i]);
4128 elementResidual_p[i] += H_s*
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
4129 elementResidual_p[i] -= H_s*
DM*
ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),p_test_dS[i]);
4130 globalConservationError += H_s*
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
4132 for (
int i=0;i<nDOF_v_test_element;i++)
4134 elementResidual_u[i] += H_s*(
ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
4135 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
4136 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
4137 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
4138 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4143 sdInfo_u_u_rowptr.data(),
4144 sdInfo_u_u_colind.data(),
4145 mom_uu_diff_ten_ext,
4146 &vel_grad_test_dS[i*nSpace])+
4147 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
4148 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4153 sdInfo_u_v_rowptr.data(),
4154 sdInfo_u_v_colind.data(),
4155 mom_uv_diff_ten_ext,
4156 &vel_grad_test_dS[i*nSpace]));
4157 elementResidual_v[i] += H_s*(
ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
4158 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
4159 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
4160 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
4161 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4166 sdInfo_v_u_rowptr.data(),
4167 sdInfo_v_u_colind.data(),
4168 mom_vu_diff_ten_ext,
4169 &vel_grad_test_dS[i*nSpace])+
4170 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
4171 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4176 sdInfo_v_v_rowptr.data(),
4177 sdInfo_v_v_colind.data(),
4178 mom_vv_diff_ten_ext,
4179 &vel_grad_test_dS[i*nSpace]));
4186 for (
int i=0;i<nDOF_test_element;i++)
4188 int eN_i = eN*nDOF_test_element+i;
4190 elementResidual_p_save.data()[eN_i] += elementResidual_p[i];
4191 mesh_volume_conservation_weak += elementResidual_mesh[i];
4192 globalResidual.data()[offset_p+stride_p*rp_l2g.data()[eN_i]]+=elementResidual_p[i];
4194 for (
int i=0;i<nDOF_v_test_element;i++)
4196 int eN_i = eN*nDOF_v_test_element+i;
4197 globalResidual.data()[offset_u+stride_u*rvel_l2g.data()[eN_i]]+=elementResidual_u[i];
4198 globalResidual.data()[offset_v+stride_v*rvel_l2g.data()[eN_i]]+=elementResidual_v[i];
4201 if (normalize_pressure)
4203 double send[4]={pa_dv,p_dv,total_volume, total_surface_area}, recv[4]={0.,0.,0.,0.};
4204 MPI_Allreduce(send, recv,4,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
4207 total_volume = recv[2];
4208 total_surface_area = recv[3];
4221 int nDOF_pressure=0;
4222 for(
int eN=0;eN<nElements_global;eN++)
4224 for (
int i=0;i<nDOF_test_element;i++)
4226 int eN_i = eN*nDOF_test_element+i;
4227 if (p_l2g.data()[eN_i] > nDOF_pressure)
4228 nDOF_pressure=p_l2g.data()[eN_i];
4232 assert(p_dof.shape(0) == nDOF_pressure);
4234 for (
int I=0;I<nDOF_pressure;I++)
4235 p_dof.data()[I] += (pa_dv - p_dv)/total_volume;
4236 double p_dv_new=0.0, pa_dv_new=0.0;
4240 for (
int eN=0 ; eN < nElements_owned ; ++eN)
4242 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
4243 for (
int j=0;j<nDOF_mesh_trial_element;j++)
4245 int eN_j = eN*nDOF_mesh_trial_element+j;
4246 element_phi[j] = phi_nodes.data()[p_l2g.data()[eN_j]];
4247 element_phi_s[j] = phi_solid_nodes.data()[p_l2g.data()[eN_j]];
4249 double element_nodes[nDOF_mesh_trial_element*3];
4250 for (
int i=0;i<nDOF_mesh_trial_element;i++)
4252 int eN_i=eN*nDOF_mesh_trial_element+i;
4253 for(
int I=0;I<3;I++)
4254 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
4256 int icase_s =
gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
4257 for (
int k=0 ; k < nQuadraturePoints_element ; ++k)
4259 int eN_k = eN*nQuadraturePoints_element + k;
4260 int eN_nDOF_trial_element = eN*nDOF_trial_element;
4262 double jac[nSpace*nSpace];
4263 double jacInv[nSpace*nSpace];
4264 double p=0.0,
pe=0.0;
4265 double jacDet, x, y,
z, dV, h_phi;
4267 double H_s =
gf_s.
H(0.,0.);
4268 ck.calculateMapping_element(eN,
4272 mesh_trial_ref.data(),
4273 mesh_grad_trial_ref.data(),
4278 dV = fabs(jacDet)*dV_ref.data()[k];
4279 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p);
4280 if (isActiveElement[eN])
4282 p_dv_new += p*H_s*dV;
4283 pa_dv_new += q_u_0.data()[eN_k]*H_s*dV;
4284 pe = p-q_u_0.data()[eN_k];
4285 p_L1 += fabs(
pe)*H_s*dV;
4286 p_L2 +=
pe*
pe*H_s*dV;
4287 if (fabs(
pe) > p_LI)
4293 assert(errors.shape(0)*errors.shape(1) == 15);
4294 MPI_Allreduce(MPI_IN_PLACE, errors.data(),(errors.shape(0)-1)*errors.shape(1),MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
4295 MPI_Allreduce(MPI_IN_PLACE, errors.data()+(errors.shape(0)-1)*errors.shape(1),1*errors.shape(1),MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
4296 assert(p_L2 >= 0.0);
4297 assert(u_L2 >= 0.0);
4298 assert(v_L2 >= 0.0);
4299 assert(velocity_L2 >= 0.0);
4303 velocity_L2 = sqrt(velocity_L2);
4308 double NONCONSERVATIVE_FORM = args.
scalar<
double>(
"NONCONSERVATIVE_FORM");
4309 double MOMENTUM_SGE = args.
scalar<
double>(
"MOMENTUM_SGE");
4310 double PRESSURE_SGE = args.
scalar<
double>(
"PRESSURE_SGE");
4311 double VELOCITY_SGE = args.
scalar<
double>(
"VELOCITY_SGE");
4312 double PRESSURE_PROJECTION_STABILIZATION = args.
scalar<
double>(
"PRESSURE_PROJECTION_STABILIZATION");
4313 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4314 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4315 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4316 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
4317 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
4318 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4319 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
4320 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4321 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4322 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4323 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4324 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4325 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4326 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4327 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4328 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4329 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4330 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4331 xt::pyarray<double>& xb_ref = args.
array<
double>(
"xb_ref");
4332 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4333 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4334 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4335 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4336 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4337 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4338 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4339 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4340 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4341 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4342 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4343 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
4344 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4345 xt::pyarray<double>& elementBoundaryDiameter = args.
array<
double>(
"elementBoundaryDiameter");
4346 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4347 double hFactor = args.
scalar<
double>(
"hFactor");
4348 int nElements_global = args.
scalar<
int>(
"nElements_global");
4349 double useRBLES = args.
scalar<
double>(
"useRBLES");
4350 double useMetrics = args.
scalar<
double>(
"useMetrics");
4351 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4352 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4353 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4354 double sigma = args.
scalar<
double>(
"sigma");
4359 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4360 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4361 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4362 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4363 double C_dg = args.
scalar<
double>(
"C_dg");
4364 double C_b = args.
scalar<
double>(
"C_b");
4365 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
4366 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
4367 const xt::pyarray<double>& eps_porous = args.
array<
double>(
"eps_porous");
4368 const xt::pyarray<double>& phi_porous = args.
array<
double>(
"phi_porous");
4369 const xt::pyarray<double>& q_velocity_porous = args.
array<
double>(
"q_velocity_porous");
4370 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
4371 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
4372 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
4373 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
4374 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
4375 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
4376 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
4377 const double LAG_LES = args.
scalar<
double>(
"LAG_LES");
4378 xt::pyarray<double> & q_eddy_viscosity_last = args.
array<
double>(
"q_eddy_viscosity_last");
4379 xt::pyarray<double> & ebqe_eddy_viscosity_last = args.
array<
double>(
"ebqe_eddy_viscosity_last");
4380 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4381 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4382 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4383 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4384 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4385 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4386 xt::pyarray<double>& p_old_dof = args.
array<
double>(
"p_old_dof");
4387 xt::pyarray<double>& u_old_dof = args.
array<
double>(
"u_old_dof");
4388 xt::pyarray<double>& v_old_dof = args.
array<
double>(
"v_old_dof");
4389 xt::pyarray<double>& w_old_dof = args.
array<
double>(
"w_old_dof");
4390 xt::pyarray<double>& g = args.
array<
double>(
"g");
4391 const double useVF = args.
scalar<
double>(
"useVF");
4392 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4393 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4394 xt::pyarray<double>& phi_nodes = args.
array<
double>(
"phi_nodes");
4395 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4396 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4397 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4398 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4399 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4400 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
4401 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
4402 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4403 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4404 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4405 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4406 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4407 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4408 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4409 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4410 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4411 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4412 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4413 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4414 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4415 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4416 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4417 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4418 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4419 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4420 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4421 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4422 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4423 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4424 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4425 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
4426 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
4427 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
4428 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
4429 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
4430 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
4431 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
4432 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
4433 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
4434 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
4435 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
4436 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
4437 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
4438 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
4439 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
4440 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
4441 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
4442 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
4443 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
4444 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
4445 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
4446 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
4447 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
4448 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
4449 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
4450 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
4451 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
4452 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
4453 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
4454 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
4455 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
4456 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
4457 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
4458 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4459 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4460 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4461 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4462 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4463 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
4464 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4465 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
4466 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4467 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4468 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
4469 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
4470 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
4471 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4472 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4473 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4474 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4475 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4476 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4477 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4478 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4479 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4480 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4481 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4482 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4483 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4484 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4485 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4486 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4487 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4488 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4489 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4490 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4491 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4492 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4493 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4494 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
4495 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
4496 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
4497 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
4498 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
4499 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
4500 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
4501 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
4502 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
4503 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
4504 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
4505 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
4506 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
4507 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
4508 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
4509 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
4510 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
4511 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
4512 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
4513 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
4514 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
4515 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
4516 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
4517 xt::pyarray<double>& ball_density = args.
array<
double>(
"ball_density");
4518 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
4519 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
4520 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
4521 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
4522 xt::pyarray<double>& ebqe_phi_s = args.
array<
double>(
"ebqe_phi_s");
4523 xt::pyarray<double>& ebq_global_grad_phi_s = args.
array<
double>(
"ebq_global_grad_phi_s");
4524 xt::pyarray<double>& ebq_particle_velocity_s = args.
array<
double>(
"ebq_particle_velocity_s");
4525 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
4526 xt::pyarray<double>& distance_to_solids = args.
array<
double>(
"distance_to_solids");
4527 xt::pyarray<int>& isActiveElement = args.
array<
int>(
"isActiveElement");
4528 xt::pyarray<int>& isActiveElement_last = args.
array<
int>(
"isActiveElement_last");
4529 int nParticles = args.
scalar<
int>(
"nParticles");
4530 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
4531 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
4532 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
4533 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
4534 double particle_beta = args.
scalar<
double>(
"particle_beta");
4535 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
4536 double ghost_penalty_constant = args.
scalar<
double>(
"ghost_penalty_constant");
4537 const bool useExact = args.
scalar<
int>(
"useExact");
4538 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
4539 std::valarray<double> particle_surfaceArea_tmp(nParticles), particle_surfaceArea_projected_tmp(nParticles), projection_direction_tmp(2), particle_volume_tmp(nParticles), particle_netForces_tmp(nParticles*3*3), particle_netMoments_tmp(nParticles*3);
4546 for(
int eN=0;eN<nElements_global;eN++)
4548 int particle_index=0;
4549 double eps_rho,eps_mu;
4551 double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
4552 elementJacobian_p_u[nDOF_test_element][nDOF_v_trial_element],
4553 elementJacobian_p_v[nDOF_test_element][nDOF_v_trial_element],
4554 elementJacobian_p_w[nDOF_test_element][nDOF_v_trial_element],
4555 elementJacobian_u_p[nDOF_v_test_element][nDOF_trial_element],
4556 elementJacobian_u_u[nDOF_v_test_element][nDOF_v_trial_element],
4557 elementJacobian_u_v[nDOF_v_test_element][nDOF_v_trial_element],
4558 elementJacobian_u_w[nDOF_v_test_element][nDOF_v_trial_element],
4559 elementJacobian_v_p[nDOF_v_test_element][nDOF_trial_element],
4560 elementJacobian_v_u[nDOF_v_test_element][nDOF_v_trial_element],
4561 elementJacobian_v_v[nDOF_v_test_element][nDOF_v_trial_element],
4562 elementJacobian_v_w[nDOF_v_test_element][nDOF_v_trial_element],
4563 elementJacobian_w_p[nDOF_v_test_element][nDOF_trial_element],
4564 elementJacobian_w_u[nDOF_v_test_element][nDOF_v_trial_element],
4565 elementJacobian_w_v[nDOF_v_test_element][nDOF_v_trial_element],
4566 elementJacobian_w_w[nDOF_v_test_element][nDOF_v_trial_element];
4567 for (
int i=0;i<nDOF_test_element;i++)
4568 for (
int j=0;j<nDOF_trial_element;j++)
4570 elementJacobian_p_p[i][j]=0.0;
4572 for (
int i=0;i<nDOF_test_element;i++)
4573 for (
int j=0;j<nDOF_v_trial_element;j++)
4575 elementJacobian_p_u[i][j]=0.0;
4576 elementJacobian_p_v[i][j]=0.0;
4577 elementJacobian_p_w[i][j]=0.0;
4578 elementJacobian_u_p[j][i]=0.0;
4579 elementJacobian_v_p[j][i]=0.0;
4580 elementJacobian_w_p[j][i]=0.0;
4582 for (
int i=0;i<nDOF_v_test_element;i++)
4583 for (
int j=0;j<nDOF_v_trial_element;j++)
4585 elementJacobian_u_u[i][j]=0.0;
4586 elementJacobian_u_v[i][j]=0.0;
4587 elementJacobian_u_w[i][j]=0.0;
4588 elementJacobian_v_u[i][j]=0.0;
4589 elementJacobian_v_v[i][j]=0.0;
4590 elementJacobian_v_w[i][j]=0.0;
4591 elementJacobian_w_u[i][j]=0.0;
4592 elementJacobian_w_v[i][j]=0.0;
4593 elementJacobian_w_w[i][j]=0.0;
4595 if(use_ball_as_particle==1 && nParticles > 0)
4597 double min_d = 1e10;
4599 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4602 mesh_dof.data()[3*mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]+0],
4603 mesh_dof.data()[3*mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]+1],
4604 mesh_dof.data()[3*mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]+2],
4605 phi_solid_nodes.data()[mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]]);
4606 if (phi_solid_nodes.data()[mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]] < min_d)
4608 min_d = phi_solid_nodes.data()[mesh_l2g.data()[eN*nDOF_mesh_trial_element+I]];
4609 particle_index = index;
4617 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
4618 for (
int j=0;j<nDOF_mesh_trial_element;j++)
4620 int eN_j = eN*nDOF_mesh_trial_element+j;
4621 element_phi[j] = phi_nodes.data()[p_l2g.data()[eN_j]];
4622 element_phi_s[j] = phi_solid_nodes.data()[p_l2g.data()[eN_j]];
4624 double element_nodes[nDOF_mesh_trial_element*3];
4625 for (
int i=0;i<nDOF_mesh_trial_element;i++)
4627 int eN_i=eN*nDOF_mesh_trial_element+i;
4628 for(
int I=0;I<3;I++)
4629 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
4631 int icase_s =
gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
4633 int icase_p =
gf_p.
calculate(element_phi, element_nodes, x_ref.data(), -
rho_1*g.data()[1], -
rho_0*g.data()[1],
false,
true);
4636 int icase_p =
gf_p.
calculate(element_phi, element_nodes, x_ref.data(), 1.,1.,
false,
false);
4637 int icase =
gf.
calculate(element_phi, element_nodes, x_ref.data(), 1.,1.,
false,
false);
4639 for (
int fluid_phase=0;fluid_phase < 2 - abs(icase); fluid_phase++)
4641 for (
int k=0;k<nQuadraturePoints_element;k++)
4643 int eN_k = eN*nQuadraturePoints_element+k,
4644 eN_k_nSpace = eN_k*nSpace,
4646 eN_nDOF_trial_element = eN*nDOF_trial_element,
4647 eN_nDOF_v_trial_element = eN*nDOF_v_trial_element;
4650 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
4652 p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0,
4680 mom_uu_diff_ten[nSpace]=
ZEROVEC,
4681 mom_vv_diff_ten[nSpace]=
ZEROVEC,
4682 mom_ww_diff_ten[nSpace]=
ZEROVEC,
4693 dmom_u_ham_grad_p[nSpace]=
ZEROVEC,
4694 dmom_u_ham_grad_u[nSpace]=
ZEROVEC,
4695 dmom_u_ham_grad_v[nSpace]=
ZEROVEC,
4700 dmom_v_ham_grad_p[nSpace]=
ZEROVEC,
4701 dmom_v_ham_grad_u[nSpace]=
ZEROVEC,
4702 dmom_v_ham_grad_v[nSpace]=
ZEROVEC,
4707 dmom_w_ham_grad_p[nSpace]=
ZEROVEC,
4708 dmom_w_ham_grad_w[nSpace]=
ZEROVEC,
4722 dpdeResidual_p_u[nDOF_v_trial_element],dpdeResidual_p_v[nDOF_v_trial_element],dpdeResidual_p_w[nDOF_v_trial_element],
4723 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_v_trial_element],
4724 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_v_trial_element],
4725 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_v_trial_element],
4726 Lstar_u_p[nDOF_test_element],
4727 Lstar_v_p[nDOF_test_element],
4728 Lstar_w_p[nDOF_test_element],
4729 Lstar_u_u[nDOF_v_test_element],
4730 Lstar_v_v[nDOF_v_test_element],
4731 Lstar_w_w[nDOF_v_test_element],
4732 Lstar_p_u[nDOF_v_test_element],
4733 Lstar_p_v[nDOF_v_test_element],
4734 Lstar_p_w[nDOF_v_test_element],
4739 dsubgridError_p_u[nDOF_v_trial_element],
4740 dsubgridError_p_v[nDOF_v_trial_element],
4741 dsubgridError_p_w[nDOF_v_trial_element],
4742 dsubgridError_u_p[nDOF_trial_element],
4743 dsubgridError_u_u[nDOF_v_trial_element],
4744 dsubgridError_v_p[nDOF_trial_element],
4745 dsubgridError_v_v[nDOF_v_trial_element],
4746 dsubgridError_w_p[nDOF_trial_element],
4747 dsubgridError_w_w[nDOF_v_trial_element],
4748 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
4749 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
4752 jacInv[nSpace*nSpace],
4753 p_trial[nDOF_trial_element], vel_trial[nDOF_v_trial_element],
4754 p_grad_trial_ib[nDOF_trial_element*nSpace], vel_grad_trial_ib[nDOF_v_trial_element*nSpace],
4755 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_v_trial_element*nSpace],
4757 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_v_test_element],
4758 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_v_test_element*nSpace],
4763 dmom_u_source[nSpace]=
ZEROVEC,
4764 dmom_v_source[nSpace]=
ZEROVEC,
4765 dmom_w_source[nSpace]=
ZEROVEC,
4768 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace]=
ZEROVEC, dmom_adv_sge[nSpace]=
ZEROVEC, dmom_ham_grad_sge[nSpace]=
ZEROVEC,
4774 dmom_u_source_s[nSpace]=
ZEROVEC,
4775 dmom_v_source_s[nSpace]=
ZEROVEC,
4776 dmom_w_source_s[nSpace]=
ZEROVEC,
4780 dmom_u_adv_u_s[nSpace]=
ZEROVEC,
4781 dmom_v_adv_v_s[nSpace]=
ZEROVEC,
4782 dmom_w_adv_w_s[nSpace]=
ZEROVEC,
4784 dmom_u_ham_grad_u_s[nSpace]=
ZEROVEC,
4785 dmom_u_ham_grad_v_s[nSpace]=
ZEROVEC,
4790 dmom_v_ham_grad_u_s[nSpace]=
ZEROVEC,
4791 dmom_v_ham_grad_v_s[nSpace]=
ZEROVEC,
4796 dmom_w_ham_grad_w_s[nSpace]=
ZEROVEC,
4808 ck.calculateMapping_element(eN,
4812 mesh_trial_ref.data(),
4813 mesh_grad_trial_ref.data(),
4818 ck.calculateH_element(eN,
4820 nodeDiametersArray.data(),
4822 mesh_trial_ref.data(),
4824 ck.calculateMappingVelocity_element(eN,
4826 mesh_velocity_dof.data(),
4828 mesh_trial_ref.data(),
4833 dV = fabs(jacDet)*dV_ref.data()[k];
4834 ck.calculateG(jacInv,G,G_dd_G,tr_G);
4837 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
4838 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
4840 ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
4841 ck_v.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_v_trial_element*nSpace],jacInv,vel_grad_trial);
4842 for (
int i=0; i < nDOF_trial_element; i++)
4844 p_trial[i] = p_trial_ref.data()[k*nDOF_trial_element + i];
4845 p_grad_trial_ib[i*nSpace + 0] = p_grad_trial[i*nSpace+0];
4846 p_grad_trial_ib[i*nSpace + 1] = p_grad_trial[i*nSpace+1];
4848 for (
int i=0; i < nDOF_v_trial_element; i++)
4850 vel_trial[i] = vel_trial_ref.data()[k*nDOF_v_trial_element + i];
4851 vel_grad_trial_ib[i*nSpace + 0] = vel_grad_trial[i*nSpace+0];
4852 vel_grad_trial_ib[i*nSpace + 1] = vel_grad_trial[i*nSpace+1];
4857 for (
int i=0; i < nDOF_trial_element; i++)
4859 if (fluid_phase == 0)
4861 if (not std::isnan(
gf_p.
VA(i)))
4864 p_grad_trial_ib[i*nSpace + 0] =
gf_p.
VA_x(i);
4865 p_grad_trial_ib[i*nSpace + 1] =
gf_p.
VA_y(i);
4870 if (not std::isnan(
gf_p.
VB(i)))
4873 p_grad_trial_ib[i*nSpace + 0] =
gf_p.
VB_x(i);
4874 p_grad_trial_ib[i*nSpace + 1] =
gf_p.
VB_y(i);
4878 if(nDOF_v_trial_element == nDOF_trial_element)
4880 for (
int vi=0; vi < nDOF_v_trial_element; vi++)
4882 if (fluid_phase == 0)
4884 if (not std::isnan(
gf.
VA(vi)))
4886 vel_trial[vi] =
gf.
VA(vi);
4887 vel_grad_trial_ib[vi*nSpace + 0] =
gf.
VA_x(vi);
4888 vel_grad_trial_ib[vi*nSpace + 1] =
gf.
VA_y(vi);
4893 if (not std::isnan(
gf.
VB(vi)))
4895 vel_trial[vi] =
gf.
VB(vi);
4896 vel_grad_trial_ib[vi*nSpace + 0] =
gf.
VB_x(vi);
4897 vel_grad_trial_ib[vi*nSpace + 1] =
gf.
VB_y(vi);
4905 for (
int vi=0; vi < nDOF_v_trial_element; vi++)
4908 if (fabs(p_trial_ref.data()[k*nDOF_trial_element + vi] - p_trial[vi]) > 1.0e-8)
4910 for (
int vj=0; vj < nDOF_trial_element; vj++)
4911 std::cout<<
"Trial "<<p_trial_ref.data()[k*nDOF_trial_element + vj]<<
'\t'<<
gf_p.
VA(vj)<<
'\t'<<
gf_p.
VB(vj)<<std::endl;
4914 if (fabs(p_grad_trial[vi*nSpace + 0] - p_grad_trial_ib[vi*nSpace+0]) > 1.0e-8)
4916 for (
int vj=0; vj < nDOF_trial_element; vj++)
4917 std::cout<<
"Grad Trial x"<<p_grad_trial[vj*nSpace + 0]<<
'\t'<<
gf_p.
VA_x(vj)<<
'\t'<<
gf_p.
VB_x(vj)<<std::endl;
4920 if (fabs(p_grad_trial[vi*nSpace + 1] - p_grad_trial_ib[vi*nSpace+1]) > 1.0e-8)
4922 for (
int vj=0; vj < nDOF_trial_element; vj++)
4923 std::cout<<
"Grad Trial y "<<p_grad_trial[vj*nSpace + 1]<<
'\t'<<
gf_p.
VA_y(vj)<<
'\t'<<
gf_p.
VB_y(vj)<<std::endl;
4927 if (fabs(vel_trial_ref.data()[k*nDOF_v_trial_element + vi] - vel_trial[vi]) > 1.0e-8)
4929 for (
int vj=0; vj < nDOF_v_trial_element; vj++)
4930 std::cout<<
"Trial "<<vel_trial_ref.data()[k*nDOF_v_trial_element + vj]<<
'\t'<<
gf.
VA(vj)<<
'\t'<<
gf.
VB(vj)<<std::endl;
4933 if (fabs(vel_grad_trial[vi*nSpace + 0] - vel_grad_trial_ib[vi*nSpace+0]) > 1.0e-8)
4935 for (
int vj=0; vj < nDOF_v_trial_element; vj++)
4936 std::cout<<
"Grad Trial x"<<vel_grad_trial[vj*nSpace + 0]<<
'\t'<<
gf.
VA_x(vj)<<
'\t'<<
gf.
VB_x(vj)<<std::endl;
4939 if (fabs(vel_grad_trial[vi*nSpace + 1] - vel_grad_trial_ib[vi*nSpace+1]) > 1.0e-8)
4941 for (
int vj=0; vj < nDOF_v_trial_element; vj++)
4942 std::cout<<
"Grad Trial y "<<vel_grad_trial[vj*nSpace + 1]<<
'\t'<<
gf.
VA_y(vj)<<
'\t'<<
gf.
VB_y(vj)<<std::endl;
4952 ck.valFromDOF(p_dof.data(),&p_l2g[eN_nDOF_trial_element],p_trial,p);
4953 ck_v.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_trial,
u);
4954 ck_v.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_trial,
v);
4955 ck.valFromDOF(p_old_dof.data(),&p_l2g[eN_nDOF_trial_element],p_trial,p_old);
4956 ck_v.valFromDOF(u_old_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_trial,u_old);
4957 ck_v.valFromDOF(v_old_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_trial,v_old);
4959 ck.gradFromDOF(p_dof.data(),&p_l2g[eN_nDOF_trial_element],p_grad_trial_ib,grad_p);
4960 ck_v.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_u);
4961 ck_v.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_v);
4962 ck.gradFromDOF(p_dof.data(),&p_l2g[eN_nDOF_trial_element],p_grad_trial_ib,grad_p_old);
4963 ck_v.gradFromDOF(u_old_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_u_old);
4964 ck_v.gradFromDOF(v_old_dof.data(),&vel_l2g[eN_nDOF_v_trial_element],vel_grad_trial_ib,grad_v_old);
4967 for (
int j=0;j<nDOF_test_element;j++)
4969 p_test_dV[j] = p_trial[j]*dV;
4970 for (
int I=0;I<nSpace;I++)
4972 p_grad_test_dV[j*nSpace+I] = p_grad_trial_ib[j*nSpace+I]*dV;
4975 for (
int j=0;j<nDOF_v_test_element;j++)
4977 vel_test_dV[j] = vel_trial[j]*dV;
4978 for (
int I=0;I<nSpace;I++)
4980 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial_ib[j*nSpace+I]*dV;
4984 for (
int j=0;j<nDOF_test_element;j++)
4986 p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
4987 for (
int I=0;I<nSpace;I++)
4989 p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;
4992 for (
int j=0;j<nDOF_v_test_element;j++)
4994 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_v_trial_element+j]*dV;
4995 for (
int I=0;I<nSpace;I++)
4997 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
5002 double div_mesh_velocity=0.0;
5003 for (
int j=0;j<nDOF_trial_element;j++)
5005 int eN_j=eN*nDOF_trial_element+j;
5006 div_mesh_velocity +=
5007 mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+0]*p_grad_trial[j*nSpace+0] +
5008 mesh_velocity_dof.data()[mesh_l2g.data()[eN_j]*3+1]*p_grad_trial[j*nSpace+1];
5010 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last.data()[eN_k])/dV;
5013 porosity = q_porosity.data()[eN_k];
5015 double ball_n[nSpace];
5016 if (use_ball_as_particle == 1 && nParticles > 0)
5018 int ball_index=
get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x,y,
z,distance_to_solids.data()[eN_k]);
5019 get_normal_to_ith_ball(nParticles, ball_center.data(), ball_radius.data(),ball_index,x,y,
z,ball_n[0],ball_n[1]);
5028 double eddy_viscosity(0.);
5029 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5030 const double H_s =
gf_s.
H(particle_eps, phi_solid.data()[eN_k]);
5036 if (fluid_phase == 0)
5047 else if (icase == -1)
5052 else if (icase == 1)
5062 double H = (1.0-useVF)*
gf.
H(eps_rho,
phi[eN_k]) + useVF*fmin(1.0,fmax(0.0,vf[eN_k]));
5063 double ImH = (1.0-useVF)*
gf.
ImH(eps_rho,
phi[eN_k]) + useVF*(1.0-fmin(1.0,fmax(0.0,vf[eN_k])));
5072 elementDiameter.data()[eN],
5073 smagorinskyConstant,
5074 turbulenceClosureModel,
5079 &normal_phi.data()[eN_k_nSpace],
5080 kappa_phi.data()[eN_k],
5084 phi_solid.data()[eN_k],
5103 q_eddy_viscosity_last.data()[eN_k],
5159 mass_source = q_mass_source.data()[eN_k];
5165 q_dragAlpha.data()[eN_k],
5166 q_dragBeta.data()[eN_k],
5179 q_velocity_sge.data()[eN_k_nSpace+0],
5180 q_velocity_sge.data()[eN_k_nSpace+1],
5181 q_velocity_sge.data()[eN_k_nSpace+1],
5182 eps_porous.data()[elementFlags.data()[eN]],
5183 phi_porous.data()[eN_k],
5184 q_velocity_porous.data()[eN_k_nSpace+0],
5185 q_velocity_porous.data()[eN_k_nSpace+1],
5186 q_velocity_porous.data()[eN_k_nSpace+1],
5195 if (turbulenceClosureModel >= 3)
5197 const double c_mu = 0.09;
5199 turbulenceClosureModel,
5211 q_turb_var_0.data()[eN_k],
5212 q_turb_var_1.data()[eN_k],
5213 &q_turb_var_grad_0.data()[eN_k_nSpace],
5233 if (NONCONSERVATIVE_FORM > 0.0)
5235 mom_u_ham -= MOVING_DOMAIN*dmom_u_acc_u*(grad_u[0]*
xt + grad_u[1]*yt);
5236 dmom_u_ham_grad_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
5237 dmom_u_ham_grad_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
5241 mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*
xt;
5242 mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt;
5243 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
5244 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
5247 if (NONCONSERVATIVE_FORM > 0.0)
5249 mom_v_ham -= MOVING_DOMAIN*dmom_v_acc_v*(grad_v[0]*
xt + grad_v[1]*yt);
5250 dmom_v_ham_grad_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
5251 dmom_v_ham_grad_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
5255 mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*
xt;
5256 mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt;
5257 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
5258 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
5264 q_mom_u_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
5270 q_mom_v_acc_beta_bdf.data()[eN_k]*q_dV_last.data()[eN_k]/dV,
5275 if (NONCONSERVATIVE_FORM > 0.0)
5277 mom_u_acc_t *= dmom_u_acc_u;
5278 mom_v_acc_t *= dmom_v_acc_v;
5283 if (NONCONSERVATIVE_FORM > 0.0)
5285 dmom_adv_sge[0] = 0.0;
5286 dmom_adv_sge[1] = 0.0;
5287 dmom_ham_grad_sge[0] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
5288 dmom_ham_grad_sge[1] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
5292 dmom_adv_sge[0] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
5293 dmom_adv_sge[1] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
5294 dmom_ham_grad_sge[0] = 0.0;
5295 dmom_ham_grad_sge[1] = 0.0;
5297 double mv_tau[nSpace]=
ZEROVEC;
5298 mv_tau[0] = dmom_adv_sge[0] + dmom_ham_grad_sge[0];
5299 mv_tau[1] = dmom_adv_sge[1] + dmom_ham_grad_sge[1];
5303 pdeResidual_p =
ck.Advection_strong(dmass_adv_u,grad_u) +
5304 ck.Advection_strong(dmass_adv_v,grad_v) +
5305 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last.data()[eN_k])/dV - div_mesh_velocity) +
5306 ck.Reaction_strong(mass_source);
5308 pdeResidual_u =
ck.Mass_strong(mom_u_acc_t) +
5309 ck.Advection_strong(dmom_adv_sge,grad_u) +
5310 ck.Hamiltonian_strong(dmom_ham_grad_sge,grad_u) +
5311 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
5312 ck.Reaction_strong(mom_u_source) -
5313 ck.Reaction_strong(dmom_u_acc_u*
u*div_mesh_velocity);
5315 pdeResidual_v =
ck.Mass_strong(mom_v_acc_t) +
5316 ck.Advection_strong(dmom_adv_sge,grad_v) +
5317 ck.Hamiltonian_strong(dmom_ham_grad_sge,grad_v) +
5318 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
5319 ck.Reaction_strong(mom_v_source) -
5320 ck.Reaction_strong(dmom_v_acc_v*
v*div_mesh_velocity);
5323 for (
int j=0;j<nDOF_v_trial_element;j++)
5325 int j_nSpace = j*nSpace;
5326 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial_ib[j_nSpace]);
5327 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial_ib[j_nSpace]);
5328 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial[j]) +
5329 ck.HamiltonianJacobian_strong(dmom_ham_grad_sge,&vel_grad_trial_ib[j_nSpace]) +
5330 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial_ib[j_nSpace]) -
5331 ck.ReactionJacobian_strong(dmom_u_acc_u*div_mesh_velocity,vel_trial[j]);
5332 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial[j]) +
5333 ck.HamiltonianJacobian_strong(dmom_ham_grad_sge,&vel_grad_trial_ib[j_nSpace]) +
5334 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial_ib[j_nSpace]) -
5335 ck.ReactionJacobian_strong(dmom_v_acc_v*div_mesh_velocity,vel_trial[j]);
5337 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial[j]);
5338 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial[j]);
5340 for (
int j=0;j<nDOF_trial_element;j++)
5342 int j_nSpace = j*nSpace;
5343 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial_ib[j_nSpace]);
5344 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial_ib[j_nSpace]);
5348 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
5350 elementDiameter.data()[eN],
5355 dmom_u_ham_grad_p[0],
5358 q_cfl.data()[eN_k]);
5365 dmom_u_ham_grad_p[0],
5368 q_cfl.data()[eN_k]);
5370 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
5371 tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0;
5405 dmom_adv_star[0] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
5406 dmom_adv_star[1] =
inertial_term*dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
5409 for (
int i=0;i<nDOF_test_element;i++)
5411 int i_nSpace = i*nSpace;
5412 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
5413 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
5416 for (
int i=0;i<nDOF_v_test_element;i++)
5418 int i_nSpace = i*nSpace;
5419 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5420 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5421 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5422 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5424 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
5425 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
5429 dmom_u_adv_u[0] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u);
5430 dmom_u_adv_u[1] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v);
5432 dmom_v_adv_v[0] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_u);
5433 dmom_v_adv_v[1] +=
inertial_term*dmom_u_acc_u*(useRBLES*subgridError_v);
5438 double level_set_normal[nSpace];
5442 double norm_exact=0.0,norm_cut=0.0;
5443 if (use_ball_as_particle)
5445 for (
int I=0;I<nSpace;I++)
5449 norm_cut += level_set_normal[I]*level_set_normal[I];
5450 norm_exact += ball_n[I]*ball_n[I];
5455 for (
int I=0;I<nSpace;I++)
5459 norm_cut += level_set_normal[I]*level_set_normal[I];
5460 norm_exact += particle_signed_distance_normals.data()[eN_k_3d+I]*particle_signed_distance_normals.data()[eN_k_3d+I];
5463 assert(std::fabs(1.0-norm_cut) < 1.0e-8);
5464 assert(std::fabs(1.0-norm_exact) < 1.0e-8);
5466 for (
int I=0;I<nSpace;I++)
5467 level_set_normal[I]*=-1.0;
5477 if (use_ball_as_particle)
5478 for (
int I=0;I<nSpace;I++)
5479 level_set_normal[I] = ball_n[I];
5481 for (
int I=0;I<nSpace;I++)
5482 level_set_normal[I] = particle_signed_distance_normals.data()[eN_k_3d+I];
5485 NONCONSERVATIVE_FORM,
5486 eN < nElements_owned,
5490 nQuadraturePoints_global,
5491 &particle_signed_distances.data()[eN_k],
5493 &particle_velocities.data()[eN_k_3d],
5494 particle_centroids.data(),
5495 use_ball_as_particle,
5498 ball_velocity.data(),
5499 ball_angular_velocity.data(),
5500 ball_density.data(),
5502 particle_penalty_constant/h_phi,
5521 q_velocity_sge.data()[eN_k_nSpace+0],
5522 q_velocity_sge.data()[eN_k_nSpace+1],
5523 q_velocity_sge.data()[eN_k_nSpace+1],
5542 dmom_u_ham_grad_u_s,
5543 dmom_u_ham_grad_v_s,
5548 dmom_v_ham_grad_u_s,
5549 dmom_v_ham_grad_v_s,
5554 dmom_w_ham_grad_w_s,
5562 &particle_netForces_tmp[0],
5563 &particle_netMoments_tmp[0],
5564 &particle_surfaceArea_tmp[0],
5565 &particle_surfaceArea_projected_tmp[0],
5566 &projection_direction_tmp[0],
5567 &particle_volume_tmp[0]);
5573 if (fluid_phase == 0)
5574 H_f =
gf.
ImH(0.,0.);
5580 assert(fluid_phase == 0);
5583 for(
int i=0;i<nDOF_test_element;i++)
5585 int i_nSpace = i*nSpace;
5586 for(
int j=0;j<nDOF_trial_element;j++)
5588 int j_nSpace = j*nSpace;
5589 if (nDOF_test_element == nDOF_v_trial_element)
5591 elementJacobian_p_p[i][j] += H_s*H_f*((1-PRESSURE_PROJECTION_STABILIZATION)*
ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_p[i]) +
5592 (1-PRESSURE_PROJECTION_STABILIZATION)*
ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_p[i]) +
5593 PRESSURE_PROJECTION_STABILIZATION*
ck.pressureProjection_weak(mom_uu_diff_ten[1], p_trial[j], 1./3., p_test_ref.data()[k*nDOF_test_element +i],dV));
5597 for(
int i=0;i<nDOF_test_element;i++)
5599 int i_nSpace = i*nSpace;
5600 for(
int j=0;j<nDOF_v_trial_element;j++)
5602 int j_nSpace = j*nSpace;
5603 elementJacobian_p_u[i][j] += H_s*H_f*(
ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial[j],&p_grad_test_dV[i_nSpace]) +
5604 ck.MassJacobian_weak(dmass_ham_u,vel_trial[j],p_test_dV[i]));
5605 elementJacobian_p_v[i][j] += H_s*H_f*(
ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial[j],&p_grad_test_dV[i_nSpace]) +
5606 ck.MassJacobian_weak(dmass_ham_v,vel_trial[j],p_test_dV[i]));
5607 if (nDOF_test_element == nDOF_v_trial_element)
5609 elementJacobian_p_u[i][j] += H_s*H_f*(1-PRESSURE_PROJECTION_STABILIZATION)*
ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]);
5610 elementJacobian_p_v[i][j] += H_s*H_f*(1-PRESSURE_PROJECTION_STABILIZATION)*
ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]);
5614 for(
int i=0;i<nDOF_v_test_element;i++)
5616 int i_nSpace = i*nSpace;
5617 for(
int j=0;j<nDOF_trial_element;j++)
5619 int j_nSpace = j*nSpace;
5620 elementJacobian_u_p[i][j] += H_s*H_f*(
ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial_ib[j_nSpace],vel_test_dV[i])+
5621 MOMENTUM_SGE*VELOCITY_SGE*
ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_u[i]));
5622 elementJacobian_v_p[i][j] += H_s*H_f*(
ck.HamiltonianJacobian_weak(dmom_v_ham_grad_p,&p_grad_trial_ib[j_nSpace],vel_test_dV[i])+
5623 MOMENTUM_SGE*VELOCITY_SGE*
ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_v[i]));
5626 for(
int i=0;i<nDOF_v_test_element;i++)
5628 int i_nSpace = i*nSpace;
5629 for(
int j=0;j<nDOF_v_trial_element;j++)
5631 int j_nSpace = j*nSpace;
5632 elementJacobian_u_u[i][j] += H_s*H_f*(
ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial[j],vel_test_dV[i]) +
5633 ck.MassJacobian_weak(dmom_u_ham_u,vel_trial[j],vel_test_dV[i]) +
5634 ck.HamiltonianJacobian_weak(dmom_u_ham_grad_u,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5635 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial[j],&vel_grad_test_dV[i_nSpace]) +
5636 ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,&vel_grad_trial_ib[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5637 ck.ReactionJacobian_weak(dmom_u_source[0]+NONCONSERVATIVE_FORM*dmom_u_acc_u*div_mesh_velocity,vel_trial[j],vel_test_dV[i]) +
5638 MOMENTUM_SGE*PRESSURE_SGE*
ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_u[i]) +
5639 MOMENTUM_SGE*VELOCITY_SGE*
ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
5640 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&vel_grad_trial_ib[j_nSpace],&vel_grad_test_dV[i_nSpace]));
5641 elementJacobian_u_v[i][j] += H_s*H_f*(
ck.HamiltonianJacobian_weak(dmom_u_ham_grad_v,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5642 ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial[j],&vel_grad_test_dV[i_nSpace]) +
5643 ck.MassJacobian_weak(dmom_u_ham_v,vel_trial[j],vel_test_dV[i]) +
5644 ck.SimpleDiffusionJacobian_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,&vel_grad_trial_ib[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5645 ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial[j],vel_test_dV[i]) +
5646 MOMENTUM_SGE*PRESSURE_SGE*
ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_u[i]));
5647 elementJacobian_v_u[i][j] += H_s*H_f*(
ck.HamiltonianJacobian_weak(dmom_v_ham_grad_u,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5648 ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial[j],&vel_grad_test_dV[i_nSpace]) +
5649 ck.MassJacobian_weak(dmom_v_ham_u,vel_trial[j],vel_test_dV[i]) +
5650 ck.SimpleDiffusionJacobian_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,&vel_grad_trial_ib[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5651 ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial[j],vel_test_dV[i]) +
5652 MOMENTUM_SGE*PRESSURE_SGE*
ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_v[i]));
5653 elementJacobian_v_v[i][j] += H_s*H_f*(
ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial[j],vel_test_dV[i]) +
5654 ck.MassJacobian_weak(dmom_v_ham_v,vel_trial[j],vel_test_dV[i]) +
5655 ck.HamiltonianJacobian_weak(dmom_v_ham_grad_v,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5656 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial[j],&vel_grad_test_dV[i_nSpace]) +
5657 ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,&vel_grad_trial_ib[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5658 ck.ReactionJacobian_weak(dmom_v_source[1]+NONCONSERVATIVE_FORM*dmom_v_acc_v*div_mesh_velocity,vel_trial[j],vel_test_dV[i]) +
5659 MOMENTUM_SGE*PRESSURE_SGE*
ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_v[i]) +
5660 MOMENTUM_SGE*VELOCITY_SGE*
ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
5661 ck.NumericalDiffusionJacobian(q_numDiff_v_last.data()[eN_k],&vel_grad_trial_ib[j_nSpace],&vel_grad_test_dV[i_nSpace]));
5666 for(
int i=0;i<nDOF_v_test_element;i++)
5668 int i_nSpace = i*nSpace;
5669 for(
int j=0;j<nDOF_v_trial_element;j++)
5671 int j_nSpace = j*nSpace;
5672 elementJacobian_u_u[i][j] += H_f*(
ck.MassJacobian_weak(dmom_u_ham_u_s,vel_trial[j],vel_test_dV[i]) +
5673 ck.HamiltonianJacobian_weak(dmom_u_ham_grad_u_s,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5674 ck.AdvectionJacobian_weak(dmom_u_adv_u_s,vel_trial[j],&vel_grad_test_dV[i_nSpace]) +
5675 ck.ReactionJacobian_weak(dmom_u_source_s[0],vel_trial[j],vel_test_dV[i]));
5677 elementJacobian_u_v[i][j] += H_f*(
ck.HamiltonianJacobian_weak(dmom_u_ham_grad_v_s,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5678 ck.ReactionJacobian_weak(dmom_u_source_s[1],vel_trial[j],vel_test_dV[i]));
5680 elementJacobian_v_u[i][j] += H_f*(
ck.HamiltonianJacobian_weak(dmom_v_ham_grad_u_s,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5681 ck.ReactionJacobian_weak(dmom_v_source_s[0],vel_trial[j],vel_test_dV[i]));
5683 elementJacobian_v_v[i][j] += H_f*(
ck.MassJacobian_weak(dmom_v_ham_v_s,vel_trial[j],vel_test_dV[i]) +
5684 ck.HamiltonianJacobian_weak(dmom_v_ham_grad_v_s,&vel_grad_trial_ib[j_nSpace],vel_test_dV[i]) +
5685 ck.AdvectionJacobian_weak(dmom_v_adv_v_s,vel_trial[j],&vel_grad_test_dV[i_nSpace]) +
5686 ck.ReactionJacobian_weak(dmom_v_source_s[1],vel_trial[j],vel_test_dV[i]));
5695 for (
int i=0;i<nDOF_test_element;i++)
5697 int eN_i = eN*nDOF_test_element+i;
5698 for (
int j=0;j<nDOF_trial_element;j++)
5700 int eN_i_j = eN_i*nDOF_trial_element+j;
5701 globalJacobian.data()[csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_p_p.data()[eN_i_j]] += elementJacobian_p_p[i][j];
5704 for (
int i=0;i<nDOF_test_element;i++)
5706 int eN_i = eN*nDOF_test_element+i;
5707 for (
int j=0;j<nDOF_v_trial_element;j++)
5709 int eN_i_j = eN_i*nDOF_v_trial_element+j;
5710 globalJacobian.data()[csrRowIndeces_p_u.data()[eN_i] + csrColumnOffsets_p_u.data()[eN_i_j]] += elementJacobian_p_u[i][j];
5711 globalJacobian.data()[csrRowIndeces_p_v.data()[eN_i] + csrColumnOffsets_p_v.data()[eN_i_j]] += elementJacobian_p_v[i][j];
5714 for (
int i=0;i<nDOF_v_test_element;i++)
5716 int eN_i = eN*nDOF_v_test_element+i;
5717 for (
int j=0;j<nDOF_trial_element;j++)
5719 int eN_i_j = eN_i*nDOF_trial_element+j;
5720 globalJacobian.data()[csrRowIndeces_u_p.data()[eN_i] + csrColumnOffsets_u_p.data()[eN_i_j]] += elementJacobian_u_p[i][j];
5721 globalJacobian.data()[csrRowIndeces_v_p.data()[eN_i] + csrColumnOffsets_v_p.data()[eN_i_j]] += elementJacobian_v_p[i][j];
5724 for (
int i=0;i<nDOF_v_test_element;i++)
5726 int eN_i = eN*nDOF_v_test_element+i;
5727 for (
int j=0;j<nDOF_v_trial_element;j++)
5729 int eN_i_j = eN_i*nDOF_v_trial_element+j;
5730 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
5731 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_u_v.data()[eN_i_j]] += elementJacobian_u_v[i][j];
5733 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_v_u.data()[eN_i_j]] += elementJacobian_v_u[i][j];
5734 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_v_v.data()[eN_i_j]] += elementJacobian_v_v[i][j];
5741 std::map<int,double> DWp_Dn_jump,DW_Dn_jump;
5742 std::map<std::pair<int, int>,
int> p_p_nz, u_u_nz, v_v_nz;
5743 double gamma_cutfem=ghost_penalty_constant,gamma_cutfem_p=ghost_penalty_constant,h_cutfem=elementBoundaryDiameter.data()[*it];
5744 int eN_nDOF_v_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_v_trial_element;
5748 for (
int i_offset=1;i_offset<nDOF_v_trial_element;i_offset++)
5751 double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
5752 v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
5753 norm_v=fmax(norm_v,sqrt(
u*
u+
v*
v));
5755 double gamma_v_dim =
rho_0*(
nu_0 + norm_v*h_cutfem + alphaBDF*h_cutfem*h_cutfem);
5756 gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
5757 if (NONCONSERVATIVE_FORM)
5758 gamma_cutfem*=gamma_v_dim;
5760 gamma_cutfem*=(gamma_v_dim/
rho_0);
5761 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5763 double Dp_Dn_jump=0.0, Du_Dn_jump=0.0, Dv_Dn_jump=0.0,dS;
5764 for (
int eN_side=0;eN_side < 2; eN_side++)
5767 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
5768 for (
int i=0;i<nDOF_test_element;i++)
5769 DWp_Dn_jump[p_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
5770 for (
int i=0;i<nDOF_v_test_element;i++)
5771 DW_Dn_jump[vel_l2g.data()[eN*nDOF_v_test_element+i]] = 0.0;
5773 for (
int eN_side=0;eN_side < 2; eN_side++)
5776 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side],
5777 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+eN_side],
5778 eN_nDOF_trial_element = eN*nDOF_trial_element,
5779 eN_nDOF_v_trial_element = eN*nDOF_v_trial_element,
5780 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5781 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5788 jac_int[nSpace*nSpace],
5790 jacInv_int[nSpace*nSpace],
5791 boundaryJac[nSpace*(nSpace-1)],
5792 metricTensor[(nSpace-1)*(nSpace-1)],
5793 metricTensorDetSqrt,
5794 p_test_dS[nDOF_test_element],vel_test_dS[nDOF_v_test_element],
5795 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_v_trial_element*nSpace],
5796 p_grad_test_dS[nDOF_trial_element*nSpace],vel_grad_test_dS[nDOF_v_trial_element*nSpace],
5797 normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
5798 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
5799 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;
5801 ck.calculateMapping_elementBoundary(eN,
5807 mesh_trial_trace_ref.data(),
5808 mesh_grad_trial_trace_ref.data(),
5809 boundaryJac_ref.data(),
5815 metricTensorDetSqrt,
5820 ck.calculateMappingVelocity_elementBoundary(eN,
5824 mesh_velocity_dof.data(),
5826 mesh_trial_trace_ref.data(),
5827 xt_int,yt_int,zt_int,
5832 dS = metricTensorDetSqrt*dS_ref.data()[kb];
5835 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,p_grad_trial_trace);
5836 ck_v.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_v_trial_element],jacInv_int,vel_grad_trial_trace);
5837 for (
int i=0;i<nDOF_test_element;i++)
5839 int eN_i = eN*nDOF_test_element + i;
5840 for (
int I=0;I<nSpace;I++)
5841 DWp_Dn_jump[p_l2g.data()[eN_i]] += p_grad_trial_trace[i*nSpace+I]*normal[I];
5843 for (
int i=0;i<nDOF_v_test_element;i++)
5845 int eN_i = eN*nDOF_v_test_element + i;
5846 for (
int I=0;I<nSpace;I++)
5847 DW_Dn_jump[vel_l2g.data()[eN_i]] += vel_grad_trial_trace[i*nSpace+I]*normal[I];
5850 for (
int eN_side=0;eN_side < 2; eN_side++)
5853 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
5854 for (
int i=0;i<nDOF_test_element;i++)
5856 int eN_i = eN*nDOF_test_element+i;
5857 for (
int eN_side2=0;eN_side2 < 2; eN_side2++)
5859 int eN2 = elementBoundaryElementsArray.data()[ebN*2+eN_side2];
5860 for (
int j=0;j<nDOF_test_element;j++)
5862 int eN_i_j = eN_i*nDOF_test_element + j;
5863 int eN2_j = eN2*nDOF_test_element + j;
5867 i*nDOF_trial_element +
5869 std::pair<int,int> ij = std::make_pair(p_l2g.data()[eN_i], p_l2g.data()[eN2_j]);
5870 if (p_p_nz.count(ij))
5872 assert(p_p_nz[ij] == csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_eb_p_p.data()[ebN_i_j]);
5875 p_p_nz[ij] = csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_eb_p_p.data()[ebN_i_j];
5879 for (
int i=0;i<nDOF_v_test_element;i++)
5881 int eN_i = eN*nDOF_v_test_element+i;
5882 for (
int eN_side2=0;eN_side2 < 2; eN_side2++)
5884 int eN2 = elementBoundaryElementsArray.data()[ebN*2+eN_side2];
5885 for (
int j=0;j<nDOF_v_test_element;j++)
5887 int eN_i_j = eN_i*nDOF_v_test_element + j;
5888 int eN2_j = eN2*nDOF_v_test_element + j;
5892 i*nDOF_v_trial_element +
5894 std::pair<int,int> ij = std::make_pair(vel_l2g.data()[eN_i], vel_l2g.data()[eN2_j]);
5895 if (u_u_nz.count(ij))
5897 assert(u_u_nz[ij] == csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]);
5900 u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
5901 if (v_v_nz.count(ij))
5903 assert(v_v_nz[ij] == csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j]);
5906 v_v_nz[ij] = csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j];
5911 for (std::map<int,double>::iterator Wi_it=DWp_Dn_jump.begin(); Wi_it!=DWp_Dn_jump.end(); ++Wi_it)
5912 for (std::map<int,double>::iterator Wj_it=DWp_Dn_jump.begin(); Wj_it!=DWp_Dn_jump.end(); ++Wj_it)
5914 int i_global = Wi_it->first,
5915 j_global = Wj_it->first;
5916 double DWp_Dn_jump_i = Wi_it->second,
5917 DWp_Dn_jump_j = Wj_it->second;
5918 std::pair<int,int> ij = std::make_pair(i_global, j_global);
5919 globalJacobian.data()[p_p_nz.at(ij)] += gamma_cutfem_p*h_cutfem*DWp_Dn_jump_j*DWp_Dn_jump_i*dS;
5921 for (std::map<int,double>::iterator Wi_it=DW_Dn_jump.begin(); Wi_it!=DW_Dn_jump.end(); ++Wi_it)
5922 for (std::map<int,double>::iterator Wj_it=DW_Dn_jump.begin(); Wj_it!=DW_Dn_jump.end(); ++Wj_it)
5924 int i_global = Wi_it->first,
5925 j_global = Wj_it->first;
5926 double DW_Dn_jump_i = Wi_it->second,
5927 DW_Dn_jump_j = Wj_it->second;
5928 std::pair<int,int> ij = std::make_pair(i_global, j_global);
5929 globalJacobian.data()[u_u_nz.at(ij)] += gamma_cutfem*h_cutfem*DW_Dn_jump_j*DW_Dn_jump_i*dS;
5930 globalJacobian.data()[v_v_nz.at(ij)] += gamma_cutfem*h_cutfem*DW_Dn_jump_j*DW_Dn_jump_i*dS;
5938 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
5940 int ebN = exteriorElementBoundariesArray.data()[ebNE],
5941 eN = elementBoundaryElementsArray.data()[ebN*2+0],
5942 eN_nDOF_trial_element = eN*nDOF_trial_element,
5943 eN_nDOF_v_trial_element = eN*nDOF_v_trial_element,
5944 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
5945 if (boundaryFlags[ebN] < 1)
5947 double eps_rho,eps_mu;
5948 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
5949 for (
int j=0;j<nDOF_mesh_trial_element;j++)
5951 int eN_j = eN*nDOF_mesh_trial_element+j;
5952 element_phi[j] = phi_nodes.data()[p_l2g.data()[eN_j]];
5953 element_phi_s[j] = phi_solid_nodes.data()[p_l2g.data()[eN_j]];
5955 double element_nodes[nDOF_mesh_trial_element*3];
5956 for (
int i=0;i<nDOF_mesh_trial_element;i++)
5958 int eN_i=eN*nDOF_mesh_trial_element+i;
5959 for(
int I=0;I<3;I++)
5960 element_nodes[i*3 + I] = mesh_dof[mesh_l2g.data()[eN_i]*3 + I];
5962 double mesh_dof_ref[nDOF_mesh_trial_element*3]={0.,0.,0.,1.,0.,0.,0.,1.,0.};
5963 double xb_ref_calc[nQuadraturePoints_elementBoundary*3];
5964 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5966 double x=0.0,y=0.0,
z=0.0;
5967 for (
int j=0;j<nDOF_mesh_trial_element;j++)
5969 int ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb;
5970 int ebN_local_kb_j = ebN_local_kb*nDOF_mesh_trial_element+j;
5971 x += mesh_dof_ref[j*3+0]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
5972 y += mesh_dof_ref[j*3+1]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
5973 z += mesh_dof_ref[j*3+2]*mesh_trial_trace_ref.data()[ebN_local_kb_j];
5975 xb_ref_calc[3*kb+0] = x;
5976 xb_ref_calc[3*kb+1] = y;
5977 xb_ref_calc[3*kb+2] =
z;
5979 int icase_s =
gf_s.
calculate(element_phi_s, element_nodes, xb_ref_calc,
true);
5983 int icase =
gf.
calculate(element_phi, element_nodes, xb_ref.data(), 1.0,1.0,
true,
false);
5985 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5987 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
5988 ebNE_kb_nSpace = ebNE_kb*nSpace,
5989 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5990 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5992 double phi_s_ext=0.0,
6001 p_old=0.0,u_old=0.0,v_old=0.0,w_old=0.0,
6004 dmom_u_acc_u_ext=0.0,
6006 dmom_v_acc_v_ext=0.0,
6008 dmom_w_acc_w_ext=0.0,
6010 dmass_adv_u_ext[nSpace]=
ZEROVEC,
6011 dmass_adv_v_ext[nSpace]=
ZEROVEC,
6012 dmass_adv_w_ext[nSpace]=
ZEROVEC,
6013 mom_u_adv_ext[nSpace]=
ZEROVEC,
6014 dmom_u_adv_u_ext[nSpace]=
ZEROVEC,
6015 dmom_u_adv_v_ext[nSpace]=
ZEROVEC,
6016 dmom_u_adv_w_ext[nSpace]=
ZEROVEC,
6017 mom_v_adv_ext[nSpace]=
ZEROVEC,
6018 dmom_v_adv_u_ext[nSpace]=
ZEROVEC,
6019 dmom_v_adv_v_ext[nSpace]=
ZEROVEC,
6020 dmom_v_adv_w_ext[nSpace]=
ZEROVEC,
6021 mom_w_adv_ext[nSpace]=
ZEROVEC,
6022 dmom_w_adv_u_ext[nSpace]=
ZEROVEC,
6023 dmom_w_adv_v_ext[nSpace]=
ZEROVEC,
6024 dmom_w_adv_w_ext[nSpace]=
ZEROVEC,
6025 mom_uu_diff_ten_ext[nSpace]=
ZEROVEC,
6026 mom_vv_diff_ten_ext[nSpace]=
ZEROVEC,
6027 mom_ww_diff_ten_ext[nSpace]=
ZEROVEC,
6028 mom_uv_diff_ten_ext[1],
6029 mom_uw_diff_ten_ext[1],
6030 mom_vu_diff_ten_ext[1],
6031 mom_vw_diff_ten_ext[1],
6032 mom_wu_diff_ten_ext[1],
6033 mom_wv_diff_ten_ext[1],
6034 mom_u_source_ext=0.0,
6035 mom_v_source_ext=0.0,
6036 mom_w_source_ext=0.0,
6038 dmom_u_ham_grad_p_ext[nSpace]=
ZEROVEC,
6039 dmom_u_ham_grad_u_ext[nSpace]=
ZEROVEC,
6040 dmom_u_ham_u_ext=0.0,
6041 dmom_u_ham_v_ext=0.0,
6042 dmom_u_ham_w_ext=0.0,
6044 dmom_v_ham_grad_p_ext[nSpace]=
ZEROVEC,
6045 dmom_v_ham_grad_v_ext[nSpace]=
ZEROVEC,
6046 dmom_v_ham_u_ext=0.0,
6047 dmom_v_ham_v_ext=0.0,
6048 dmom_v_ham_w_ext=0.0,
6050 dmom_w_ham_grad_p_ext[nSpace]=
ZEROVEC,
6051 dmom_w_ham_grad_w_ext[nSpace]=
ZEROVEC,
6052 dmom_w_ham_u_ext=0.0,
6053 dmom_w_ham_v_ext=0.0,
6054 dmom_w_ham_w_ext=0.0,
6055 dmom_u_adv_p_ext[nSpace]=
ZEROVEC,
6056 dmom_v_adv_p_ext[nSpace]=
ZEROVEC,
6057 dmom_w_adv_p_ext[nSpace]=
ZEROVEC,
6058 dflux_mass_u_ext=0.0,
6059 dflux_mass_v_ext=0.0,
6060 dflux_mass_w_ext=0.0,
6061 dflux_mom_u_adv_p_ext=0.0,
6062 dflux_mom_u_adv_u_ext=0.0,
6063 dflux_mom_u_adv_v_ext=0.0,
6064 dflux_mom_u_adv_w_ext=0.0,
6065 dflux_mom_v_adv_p_ext=0.0,
6066 dflux_mom_v_adv_u_ext=0.0,
6067 dflux_mom_v_adv_v_ext=0.0,
6068 dflux_mom_v_adv_w_ext=0.0,
6069 dflux_mom_w_adv_p_ext=0.0,
6070 dflux_mom_w_adv_u_ext=0.0,
6071 dflux_mom_w_adv_v_ext=0.0,
6072 dflux_mom_w_adv_w_ext=0.0,
6077 bc_mom_u_acc_ext=0.0,
6078 bc_dmom_u_acc_u_ext=0.0,
6079 bc_mom_v_acc_ext=0.0,
6080 bc_dmom_v_acc_v_ext=0.0,
6081 bc_mom_w_acc_ext=0.0,
6082 bc_dmom_w_acc_w_ext=0.0,
6083 bc_mass_adv_ext[nSpace]=
ZEROVEC,
6084 bc_dmass_adv_u_ext[nSpace]=
ZEROVEC,
6085 bc_dmass_adv_v_ext[nSpace]=
ZEROVEC,
6086 bc_dmass_adv_w_ext[nSpace]=
ZEROVEC,
6087 bc_mom_u_adv_ext[nSpace]=
ZEROVEC,
6088 bc_dmom_u_adv_u_ext[nSpace]=
ZEROVEC,
6089 bc_dmom_u_adv_v_ext[nSpace]=
ZEROVEC,
6090 bc_dmom_u_adv_w_ext[nSpace]=
ZEROVEC,
6091 bc_mom_v_adv_ext[nSpace]=
ZEROVEC,
6092 bc_dmom_v_adv_u_ext[nSpace]=
ZEROVEC,
6093 bc_dmom_v_adv_v_ext[nSpace]=
ZEROVEC,
6094 bc_dmom_v_adv_w_ext[nSpace]=
ZEROVEC,
6095 bc_mom_w_adv_ext[nSpace]=
ZEROVEC,
6096 bc_dmom_w_adv_u_ext[nSpace]=
ZEROVEC,
6097 bc_dmom_w_adv_v_ext[nSpace]=
ZEROVEC,
6098 bc_dmom_w_adv_w_ext[nSpace]=
ZEROVEC,
6099 bc_mom_uu_diff_ten_ext[nSpace]=
ZEROVEC,
6100 bc_mom_vv_diff_ten_ext[nSpace]=
ZEROVEC,
6101 bc_mom_ww_diff_ten_ext[nSpace]=
ZEROVEC,
6102 bc_mom_uv_diff_ten_ext[1],
6103 bc_mom_uw_diff_ten_ext[1],
6104 bc_mom_vu_diff_ten_ext[1],
6105 bc_mom_vw_diff_ten_ext[1],
6106 bc_mom_wu_diff_ten_ext[1],
6107 bc_mom_wv_diff_ten_ext[1],
6108 bc_mom_u_source_ext=0.0,
6109 bc_mom_v_source_ext=0.0,
6110 bc_mom_w_source_ext=0.0,
6111 bc_mom_u_ham_ext=0.0,
6112 bc_dmom_u_ham_grad_p_ext[nSpace]=
ZEROVEC,
6113 bc_dmom_u_ham_grad_u_ext[nSpace]=
ZEROVEC,
6114 bc_dmom_u_ham_u_ext=0.0,
6115 bc_dmom_u_ham_v_ext=0.0,
6116 bc_dmom_u_ham_w_ext=0.0,
6117 bc_mom_v_ham_ext=0.0,
6118 bc_dmom_v_ham_grad_p_ext[nSpace]=
ZEROVEC,
6119 bc_dmom_v_ham_grad_v_ext[nSpace]=
ZEROVEC,
6120 bc_dmom_v_ham_u_ext=0.0,
6121 bc_dmom_v_ham_v_ext=0.0,
6122 bc_dmom_v_ham_w_ext=0.0,
6123 bc_mom_w_ham_ext=0.0,
6124 bc_dmom_w_ham_grad_p_ext[nSpace]=
ZEROVEC,
6125 bc_dmom_w_ham_grad_w_ext[nSpace]=
ZEROVEC,
6126 bc_dmom_w_ham_u_ext=0.0,
6127 bc_dmom_w_ham_v_ext=0.0,
6128 bc_dmom_w_ham_w_ext=0.0,
6129 fluxJacobian_p_p[nDOF_trial_element],
6130 fluxJacobian_p_u[nDOF_v_trial_element],
6131 fluxJacobian_p_v[nDOF_v_trial_element],
6132 fluxJacobian_p_w[nDOF_v_trial_element],
6133 fluxJacobian_u_p[nDOF_trial_element],
6134 fluxJacobian_u_u[nDOF_v_trial_element],
6135 fluxJacobian_u_v[nDOF_v_trial_element],
6136 fluxJacobian_u_w[nDOF_v_trial_element],
6137 fluxJacobian_v_p[nDOF_trial_element],
6138 fluxJacobian_v_u[nDOF_v_trial_element],
6139 fluxJacobian_v_v[nDOF_v_trial_element],
6140 fluxJacobian_v_w[nDOF_v_trial_element],
6141 fluxJacobian_w_p[nDOF_trial_element],
6142 fluxJacobian_w_u[nDOF_v_trial_element],
6143 fluxJacobian_w_v[nDOF_v_trial_element],
6144 fluxJacobian_w_w[nDOF_v_trial_element],
6145 jac_ext[nSpace*nSpace],
6147 jacInv_ext[nSpace*nSpace],
6148 boundaryJac[nSpace*(nSpace-1)],
6149 metricTensor[(nSpace-1)*(nSpace-1)],
6150 metricTensorDetSqrt,
6151 p_grad_trial_trace[nDOF_trial_element*nSpace],
6152 vel_grad_trial_trace[nDOF_v_trial_element*nSpace],
6154 p_test_dS[nDOF_test_element],
6155 vel_test_dS[nDOF_v_test_element],
6157 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
6158 vel_grad_test_dS[nDOF_v_trial_element*nSpace],
6162 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
6165 ck.calculateMapping_elementBoundary(eN,
6171 mesh_trial_trace_ref.data(),
6172 mesh_grad_trial_trace_ref.data(),
6173 boundaryJac_ref.data(),
6179 metricTensorDetSqrt,
6183 ck.calculateMappingVelocity_elementBoundary(eN,
6187 mesh_velocity_dof.data(),
6189 mesh_trial_trace_ref.data(),
6190 xt_ext,yt_ext,zt_ext,
6196 dS = metricTensorDetSqrt*dS_ref.data()[kb];
6197 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
6198 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
6200 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
6201 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
6205 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
6206 ck_v.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_v_trial_element],jacInv_ext,vel_grad_trial_trace);
6208 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext);
6209 ck_v.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],u_ext);
6210 ck_v.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],v_ext);
6211 ck.valFromDOF(p_old_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_old);
6212 ck_v.valFromDOF(u_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],u_old);
6213 ck_v.valFromDOF(v_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_v_test_element],v_old);
6214 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
6215 ck_v.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_ext);
6216 ck_v.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_ext);
6217 ck.gradFromDOF(p_old_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_old);
6218 ck_v.gradFromDOF(u_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_u_old);
6219 ck_v.gradFromDOF(v_old_dof.data(),&vel_l2g.data()[eN_nDOF_v_trial_element],vel_grad_trial_trace,grad_v_old);
6220 ck.valFromDOF(phi_solid_nodes.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],phi_s_ext);
6222 for (
int j=0;j<nDOF_test_element;j++)
6224 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
6227 for (
int j=0;j<nDOF_v_test_element;j++)
6229 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_v_test_element+j]*dS;
6230 for (
int I=0;I<nSpace;I++)
6231 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
6236 bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
6238 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*(ebqe_bc_u_ext.data()[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
6239 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*(ebqe_bc_v_ext.data()[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
6241 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
6245 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);
6246 if (use_ball_as_particle == 1 && nParticles > 0)
6248 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(),x_ext,y_ext,z_ext,ebqe_phi_s.data()[ebNE_kb]);
6251 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
6253 double H = (1.0-useVF)*
gf.
H(eps_rho,ebqe_phi_ext.data()[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,ebqe_vf_ext.data()[ebNE_kb]));
6254 double ImH = (1.0-useVF)*
gf.
ImH(eps_rho,ebqe_phi_ext.data()[ebNE_kb]) + useVF*(1.0-fmin(1.0,fmax(0.0,ebqe_vf_ext.data()[ebNE_kb])));
6262 elementDiameter.data()[eN],
6263 smagorinskyConstant,
6264 turbulenceClosureModel,
6267 ebqe_vf_ext.data()[ebNE_kb],
6268 ebqe_phi_ext.data()[ebNE_kb],
6269 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
6270 ebqe_kappa_phi_ext.data()[ebNE_kb],
6274 ebqe_phi_s.data()[ebNE_kb],
6293 ebqe_eddy_viscosity_last.data()[ebNE_kb],
6316 mom_uu_diff_ten_ext,
6317 mom_vv_diff_ten_ext,
6318 mom_ww_diff_ten_ext,
6319 mom_uv_diff_ten_ext,
6320 mom_uw_diff_ten_ext,
6321 mom_vu_diff_ten_ext,
6322 mom_vw_diff_ten_ext,
6323 mom_wu_diff_ten_ext,
6324 mom_wv_diff_ten_ext,
6329 dmom_u_ham_grad_p_ext,
6330 dmom_u_ham_grad_u_ext,
6335 dmom_v_ham_grad_p_ext,
6336 dmom_v_ham_grad_v_ext,
6341 dmom_w_ham_grad_p_ext,
6342 dmom_w_ham_grad_w_ext,
6350 H = (1.0-useVF)*
gf.
H(eps_rho,bc_ebqe_phi_ext.data()[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,bc_ebqe_vf_ext.data()[ebNE_kb]));
6351 ImH = (1.0-useVF)*
gf.
ImH(eps_rho,bc_ebqe_phi_ext.data()[ebNE_kb]) + useVF*(1.0-fmin(1.0,fmax(0.0,bc_ebqe_vf_ext.data()[ebNE_kb])));
6359 elementDiameter.data()[eN],
6360 smagorinskyConstant,
6361 turbulenceClosureModel,
6364 bc_ebqe_vf_ext.data()[ebNE_kb],
6365 bc_ebqe_phi_ext.data()[ebNE_kb],
6366 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
6367 ebqe_kappa_phi_ext.data()[ebNE_kb],
6371 ebqe_phi_s.data()[ebNE_kb],
6389 bc_eddy_viscosity_ext,
6390 ebqe_eddy_viscosity_last.data()[ebNE_kb],
6392 bc_dmom_u_acc_u_ext,
6394 bc_dmom_v_acc_v_ext,
6396 bc_dmom_w_acc_w_ext,
6402 bc_dmom_u_adv_u_ext,
6403 bc_dmom_u_adv_v_ext,
6404 bc_dmom_u_adv_w_ext,
6406 bc_dmom_v_adv_u_ext,
6407 bc_dmom_v_adv_v_ext,
6408 bc_dmom_v_adv_w_ext,
6410 bc_dmom_w_adv_u_ext,
6411 bc_dmom_w_adv_v_ext,
6412 bc_dmom_w_adv_w_ext,
6413 bc_mom_uu_diff_ten_ext,
6414 bc_mom_vv_diff_ten_ext,
6415 bc_mom_ww_diff_ten_ext,
6416 bc_mom_uv_diff_ten_ext,
6417 bc_mom_uw_diff_ten_ext,
6418 bc_mom_vu_diff_ten_ext,
6419 bc_mom_vw_diff_ten_ext,
6420 bc_mom_wu_diff_ten_ext,
6421 bc_mom_wv_diff_ten_ext,
6422 bc_mom_u_source_ext,
6423 bc_mom_v_source_ext,
6424 bc_mom_w_source_ext,
6426 bc_dmom_u_ham_grad_p_ext,
6427 bc_dmom_u_ham_grad_u_ext,
6428 bc_dmom_u_ham_u_ext,
6429 bc_dmom_u_ham_v_ext,
6430 bc_dmom_u_ham_w_ext,
6432 bc_dmom_v_ham_grad_p_ext,
6433 bc_dmom_v_ham_grad_v_ext,
6434 bc_dmom_v_ham_u_ext,
6435 bc_dmom_v_ham_v_ext,
6436 bc_dmom_v_ham_w_ext,
6438 bc_dmom_w_ham_grad_p_ext,
6439 bc_dmom_w_ham_grad_w_ext,
6440 bc_dmom_w_ham_u_ext,
6441 bc_dmom_w_ham_v_ext,
6442 bc_dmom_w_ham_w_ext,
6447 if (turbulenceClosureModel >= 3)
6449 const double turb_var_grad_0_dummy[nSpace] =
ZEROVEC;
6450 const double c_mu = 0.09;
6452 turbulenceClosureModel,
6460 ebqe_vf_ext.data()[ebNE_kb],
6461 ebqe_phi_ext.data()[ebNE_kb],
6464 ebqe_turb_var_0.data()[ebNE_kb],
6465 ebqe_turb_var_1.data()[ebNE_kb],
6466 turb_var_grad_0_dummy,
6468 mom_uu_diff_ten_ext,
6469 mom_vv_diff_ten_ext,
6470 mom_ww_diff_ten_ext,
6471 mom_uv_diff_ten_ext,
6472 mom_uw_diff_ten_ext,
6473 mom_vu_diff_ten_ext,
6474 mom_vw_diff_ten_ext,
6475 mom_wu_diff_ten_ext,
6476 mom_wv_diff_ten_ext,
6482 turbulenceClosureModel,
6490 ebqe_vf_ext.data()[ebNE_kb],
6491 ebqe_phi_ext.data()[ebNE_kb],
6494 ebqe_turb_var_0.data()[ebNE_kb],
6495 ebqe_turb_var_1.data()[ebNE_kb],
6496 turb_var_grad_0_dummy,
6497 bc_eddy_viscosity_ext,
6498 bc_mom_uu_diff_ten_ext,
6499 bc_mom_vv_diff_ten_ext,
6500 bc_mom_ww_diff_ten_ext,
6501 bc_mom_uv_diff_ten_ext,
6502 bc_mom_uw_diff_ten_ext,
6503 bc_mom_vu_diff_ten_ext,
6504 bc_mom_vw_diff_ten_ext,
6505 bc_mom_wu_diff_ten_ext,
6506 bc_mom_wv_diff_ten_ext,
6507 bc_mom_u_source_ext,
6508 bc_mom_v_source_ext,
6509 bc_mom_w_source_ext);
6514 if (NONCONSERVATIVE_FORM > 0.0)
6516 mom_u_ham_ext -= MOVING_DOMAIN*dmom_u_acc_u_ext*(grad_u_ext[0]*xt_ext +
6517 grad_u_ext[1]*yt_ext);
6518 dmom_u_ham_grad_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
6519 dmom_u_ham_grad_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
6523 mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext;
6524 mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext;
6525 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
6526 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
6529 if (NONCONSERVATIVE_FORM > 0.0)
6531 mom_v_ham_ext -= MOVING_DOMAIN*dmom_v_acc_v_ext*(grad_v_ext[0]*xt_ext +
6532 grad_v_ext[1]*yt_ext);
6533 dmom_v_ham_grad_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
6534 dmom_v_ham_grad_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
6538 mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext;
6539 mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext;
6540 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
6541 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
6545 if (NONCONSERVATIVE_FORM < 1.0)
6547 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext;
6548 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext;
6550 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext;
6551 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext;
6557 isDOFBoundary_p.data()[ebNE_kb],
6558 isDOFBoundary_u.data()[ebNE_kb],
6559 isDOFBoundary_v.data()[ebNE_kb],
6560 isDOFBoundary_w.data()[ebNE_kb],
6561 isAdvectiveFluxBoundary_p.data()[ebNE_kb],
6562 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
6563 isAdvectiveFluxBoundary_v.data()[ebNE_kb],
6564 isAdvectiveFluxBoundary_w.data()[ebNE_kb],
6565 dmom_u_ham_grad_p_ext[0],
6574 ebqe_bc_flux_mass_ext.data()[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),
6575 ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
6576 ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
6577 ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
6590 dmom_u_ham_grad_u_ext,
6605 dflux_mom_u_adv_p_ext,
6606 dflux_mom_u_adv_u_ext,
6607 dflux_mom_u_adv_v_ext,
6608 dflux_mom_u_adv_w_ext,
6609 dflux_mom_v_adv_p_ext,
6610 dflux_mom_v_adv_u_ext,
6611 dflux_mom_v_adv_v_ext,
6612 dflux_mom_v_adv_w_ext,
6613 dflux_mom_w_adv_p_ext,
6614 dflux_mom_w_adv_u_ext,
6615 dflux_mom_w_adv_v_ext,
6616 dflux_mom_w_adv_w_ext);
6620 ck.calculateGScale(G,normal,h_penalty);
6621 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
6622 if (isActiveElement[eN])
6625 for (
int j=0;j<nDOF_trial_element;j++)
6627 int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
6628 fluxJacobian_p_p[j]=0.0;
6629 fluxJacobian_u_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
6630 fluxJacobian_v_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
6632 for (
int j=0;j<nDOF_v_trial_element;j++)
6634 int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_v_trial_element+j;
6635 fluxJacobian_p_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
6636 fluxJacobian_p_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
6637 fluxJacobian_u_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
6639 ebqe_phi_ext.data()[ebNE_kb],
6640 sdInfo_u_u_rowptr.data(),
6641 sdInfo_u_u_colind.data(),
6642 isDOFBoundary_u.data()[ebNE_kb],
6643 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
6645 mom_uu_diff_ten_ext,
6646 vel_trial_trace_ref.data()[ebN_local_kb_j],
6647 &vel_grad_trial_trace[j_nSpace],
6649 fluxJacobian_u_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
6651 ebqe_phi_ext.data()[ebNE_kb],
6652 sdInfo_u_v_rowptr.data(),
6653 sdInfo_u_v_colind.data(),
6654 isDOFBoundary_v.data()[ebNE_kb],
6655 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
6657 mom_uv_diff_ten_ext,
6658 vel_trial_trace_ref.data()[ebN_local_kb_j],
6659 &vel_grad_trial_trace[j_nSpace],
6662 fluxJacobian_v_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
6664 ebqe_phi_ext.data()[ebNE_kb],
6665 sdInfo_v_u_rowptr.data(),
6666 sdInfo_v_u_colind.data(),
6667 isDOFBoundary_u.data()[ebNE_kb],
6668 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
6670 mom_vu_diff_ten_ext,
6671 vel_trial_trace_ref.data()[ebN_local_kb_j],
6672 &vel_grad_trial_trace[j_nSpace],
6674 fluxJacobian_v_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
6676 ebqe_phi_ext.data()[ebNE_kb],
6677 sdInfo_v_v_rowptr.data(),
6678 sdInfo_v_v_colind.data(),
6679 isDOFBoundary_v.data()[ebNE_kb],
6680 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
6682 mom_vv_diff_ten_ext,
6683 vel_trial_trace_ref.data()[ebN_local_kb_j],
6684 &vel_grad_trial_trace[j_nSpace],
6691 const double H_s =
gf_s.
H(particle_eps, ebqe_phi_s[ebNE_kb]);
6692 if (isActiveElement[eN])
6694 for (
int i=0;i<nDOF_test_element;i++)
6696 int eN_i = eN*nDOF_test_element+i;
6697 for (
int j=0;j<nDOF_trial_element;j++)
6699 int eN_j = eN*nDOF_trial_element+j;
6700 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;
6702 globalJacobian.data()[csrRowIndeces_p_p[eN_i] + csrColumnOffsets_eb_p_p.data()[ebN_i_j]] += H_s*fluxJacobian_p_p[j]*p_test_dS[i];
6705 for (
int i=0;i<nDOF_test_element;i++)
6707 int eN_i = eN*nDOF_test_element+i;
6708 for (
int j=0;j<nDOF_v_trial_element;j++)
6710 int eN_j = eN*nDOF_v_trial_element+j;
6712 globalJacobian.data()[csrRowIndeces_p_u.data()[eN_i] + csrColumnOffsets_eb_p_u.data()[ebN_i_j]] += H_s*fluxJacobian_p_u[j]*p_test_dS[i];
6713 globalJacobian.data()[csrRowIndeces_p_v.data()[eN_i] + csrColumnOffsets_eb_p_v.data()[ebN_i_j]] += H_s*fluxJacobian_p_v[j]*p_test_dS[i];
6716 for (
int i=0;i<nDOF_v_test_element;i++)
6718 int eN_i = eN*nDOF_v_test_element+i;
6719 for (
int j=0;j<nDOF_trial_element;j++)
6722 globalJacobian.data()[csrRowIndeces_u_p.data()[eN_i] + csrColumnOffsets_eb_u_p.data()[ebN_i_j]] += H_s*fluxJacobian_u_p[j]*vel_test_dS[i];
6723 globalJacobian.data()[csrRowIndeces_v_p.data()[eN_i] + csrColumnOffsets_eb_v_p.data()[ebN_i_j]] += H_s*fluxJacobian_v_p[j]*vel_test_dS[i];
6726 for (
int i=0;i<nDOF_v_test_element;i++)
6728 int eN_i = eN*nDOF_v_test_element+i;
6729 for (
int j=0;j<nDOF_v_trial_element;j++)
6731 int eN_j = eN*nDOF_v_trial_element+j;
6733 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] +=
6734 H_s*(fluxJacobian_u_u[j]*vel_test_dS[i]+
6735 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
6736 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
6738 vel_trial_trace_ref.data()[ebN_local_kb_j],
6740 sdInfo_u_u_rowptr.data(),
6741 sdInfo_u_u_colind.data(),
6742 mom_uu_diff_ten_ext,
6743 &vel_grad_test_dS[i*nSpace]));
6744 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_eb_u_v.data()[ebN_i_j]] +=
6745 H_s*(fluxJacobian_u_v[j]*vel_test_dS[i]+
6746 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
6747 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
6749 vel_trial_trace_ref.data()[ebN_local_kb_j],
6751 sdInfo_u_v_rowptr.data(),
6752 sdInfo_u_v_colind.data(),
6753 mom_uv_diff_ten_ext,
6754 &vel_grad_test_dS[i*nSpace]));
6755 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_eb_v_u.data()[ebN_i_j]] +=
6756 H_s*(fluxJacobian_v_u[j]*vel_test_dS[i]+
6757 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
6758 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
6760 vel_trial_trace_ref.data()[ebN_local_kb_j],
6762 sdInfo_v_u_rowptr.data(),
6763 sdInfo_v_u_colind.data(),
6764 mom_vu_diff_ten_ext,
6765 &vel_grad_test_dS[i*nSpace]));
6766 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j]] +=
6767 H_s*(fluxJacobian_v_v[j]*vel_test_dS[i]+
6768 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
6769 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
6771 vel_trial_trace_ref.data()[ebN_local_kb_j],
6773 sdInfo_v_v_rowptr.data(),
6774 sdInfo_v_v_colind.data(),
6775 mom_vv_diff_ten_ext,
6776 &vel_grad_test_dS[i*nSpace]));