1581 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1582 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1583 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1584 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1585 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1586 double PSTAB = args.
scalar<
double>(
"PSTAB");
1587 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1588 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
1589 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1590 int nDOF_per_element_pressure = args.
scalar<
int>(
"nDOF_per_element_pressure");
1591 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1592 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1593 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1594 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1595 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
1596 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
1597 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
1598 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
1599 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1600 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1601 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
1602 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1603 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1604 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1605 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1606 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1607 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1608 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1609 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1610 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1611 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1612 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1613 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1614 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1615 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1616 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1617 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1618 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1619 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1620 double hFactor = args.
scalar<
double>(
"hFactor");
1621 int nElements_global = args.
scalar<
int>(
"nElements_global");
1622 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1623 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
1624 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1625 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
1626 double useRBLES = args.
scalar<
double>(
"useRBLES");
1627 double useMetrics = args.
scalar<
double>(
"useMetrics");
1628 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1629 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1630 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1631 double sigma = args.
scalar<
double>(
"sigma");
1636 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1637 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1638 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1639 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1640 double C_dc = args.
scalar<
double>(
"C_dc");
1641 double C_b = args.
scalar<
double>(
"C_b");
1642 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1643 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
1644 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
1645 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
1646 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
1647 xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1648 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
1649 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
1650 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1651 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
1652 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
1653 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1654 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1655 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1656 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1657 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1658 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1659 xt::pyarray<double>& q_eddy_viscosity = args.
array<
double>(
"q_eddy_viscosity");
1660 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1661 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1662 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1663 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1664 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1665 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1666 xt::pyarray<double>& u_dof_old = args.
array<
double>(
"u_dof_old");
1667 xt::pyarray<double>& v_dof_old = args.
array<
double>(
"v_dof_old");
1668 xt::pyarray<double>& w_dof_old = args.
array<
double>(
"w_dof_old");
1669 xt::pyarray<double>& u_dof_old_old = args.
array<
double>(
"u_dof_old_old");
1670 xt::pyarray<double>& v_dof_old_old = args.
array<
double>(
"v_dof_old_old");
1671 xt::pyarray<double>& w_dof_old_old = args.
array<
double>(
"w_dof_old_old");
1672 xt::pyarray<double>& uStar_dof = args.
array<
double>(
"uStar_dof");
1673 xt::pyarray<double>& vStar_dof = args.
array<
double>(
"vStar_dof");
1674 xt::pyarray<double>& wStar_dof = args.
array<
double>(
"wStar_dof");
1675 xt::pyarray<double>& g = args.
array<
double>(
"g");
1676 const double useVF = args.
scalar<
double>(
"useVF");
1677 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1678 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1679 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
1680 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1681 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1682 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1683 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1684 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1685 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1686 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1687 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1688 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1689 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
1690 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
1691 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1692 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
1693 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1694 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1695 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1696 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1697 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1698 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1699 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1700 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1701 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1702 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1703 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1704 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1705 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1706 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1707 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1708 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1709 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1710 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1711 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1712 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1713 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1714 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1715 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1716 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1717 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1718 int offset_p = args.
scalar<
int>(
"offset_p");
1719 int offset_u = args.
scalar<
int>(
"offset_u");
1720 int offset_v = args.
scalar<
int>(
"offset_v");
1721 int offset_w = args.
scalar<
int>(
"offset_w");
1722 int stride_p = args.
scalar<
int>(
"stride_p");
1723 int stride_u = args.
scalar<
int>(
"stride_u");
1724 int stride_v = args.
scalar<
int>(
"stride_v");
1725 int stride_w = args.
scalar<
int>(
"stride_w");
1726 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1727 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1728 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1729 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
1730 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1731 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1732 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1733 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1734 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1735 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1736 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1737 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1738 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
1739 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1740 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1741 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1742 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1743 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1744 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1745 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1746 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1747 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1748 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1749 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1750 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1751 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1752 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1753 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1754 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1755 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1756 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1757 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1758 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1759 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1760 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1761 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1762 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1763 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1764 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1765 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1766 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1767 xt::pyarray<double>& q_grad_u = args.
array<
double>(
"q_grad_u");
1768 xt::pyarray<double>& q_grad_v = args.
array<
double>(
"q_grad_v");
1769 xt::pyarray<double>& q_grad_w = args.
array<
double>(
"q_grad_w");
1770 xt::pyarray<double>& q_divU = args.
array<
double>(
"q_divU");
1771 xt::pyarray<double>& ebqe_grad_u = args.
array<
double>(
"ebqe_grad_u");
1772 xt::pyarray<double>& ebqe_grad_v = args.
array<
double>(
"ebqe_grad_v");
1773 xt::pyarray<double>& ebqe_grad_w = args.
array<
double>(
"ebqe_grad_w");
1774 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1775 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1776 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
1777 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1778 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1779 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1780 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1781 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1782 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1783 xt::pyarray<double>& q_rho = args.
array<
double>(
"q_rho");
1784 xt::pyarray<double>& ebqe_rho = args.
array<
double>(
"ebqe_rho");
1785 xt::pyarray<double>& q_nu = args.
array<
double>(
"q_nu");
1786 xt::pyarray<double>& ebqe_nu = args.
array<
double>(
"ebqe_nu");
1787 int nParticles = args.
scalar<
int>(
"nParticles");
1788 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
1789 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
1790 double particle_beta = args.
scalar<
double>(
"particle_beta");
1791 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
1792 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
1793 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
1794 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
1795 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
1796 xt::pyarray<double>& particle_netForces = args.
array<
double>(
"particle_netForces");
1797 xt::pyarray<double>& particle_netMoments = args.
array<
double>(
"particle_netMoments");
1798 xt::pyarray<double>& particle_surfaceArea = args.
array<
double>(
"particle_surfaceArea");
1799 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
1800 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
1801 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
1802 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
1803 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
1804 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
1805 xt::pyarray<double>& phisError = args.
array<
double>(
"phisError");
1806 xt::pyarray<double>& phisErrorNodal = args.
array<
double>(
"phisErrorNodal");
1807 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
1808 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
1810 double cE = args.
scalar<
double>(
"cE");
1811 int MULTIPLY_EXTERNAL_FORCE_BY_DENSITY = args.
scalar<
int>(
"MULTIPLY_EXTERNAL_FORCE_BY_DENSITY");
1812 xt::pyarray<double>& forcex = args.
array<
double>(
"forcex");
1813 xt::pyarray<double>& forcey = args.
array<
double>(
"forcey");
1814 xt::pyarray<double>& forcez = args.
array<
double>(
"forcez");
1815 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
1816 double dt = args.
scalar<
double>(
"dt");
1817 xt::pyarray<double>& quantDOFs = args.
array<
double>(
"quantDOFs");
1818 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
1819 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
1820 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
1821 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
1822 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
1823 double order_polynomial = args.
scalar<
double>(
"order_polynomial");
1824 xt::pyarray<double>& isActiveDOF = args.
array<
double>(
"isActiveDOF");
1825 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
1826 xt::pyarray<double>& ncDrag = args.
array<
double>(
"ncDrag");
1827 xt::pyarray<double>& betaDrag = args.
array<
double>(
"betaDrag");
1828 xt::pyarray<double>& vos_vel_nodes = args.
array<
double>(
"vos_vel_nodes");
1829 xt::pyarray<double>& entropyResidualPerNode = args.
array<
double>(
"entropyResidualPerNode");
1830 xt::pyarray<double>& laggedEntropyResidualPerNode = args.
array<
double>(
"laggedEntropyResidualPerNode");
1831 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
1832 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
1833 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
1834 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
1835 int NNZ_1D = args.
scalar<
int>(
"NNZ_1D");
1836 xt::pyarray<int>& csrRowIndeces_1D = args.
array<
int>(
"csrRowIndeces_1D");
1837 xt::pyarray<int>& csrColumnOffsets_1D = args.
array<
int>(
"csrColumnOffsets_1D");
1838 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
1839 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
1840 xt::pyarray<double>& isBoundary_1D = args.
array<
double>(
"isBoundary_1D");
1841 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
1847 double cut_cell_boundary_length=0.0, p_force_x=0.0, p_force_y=0.0;
1848 double element_uStar_He[nElements_global], element_vStar_He[nElements_global];
1851 den_hi.resize(numDOFs_1D,0.0);
1861 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1867 if (
psi.size() != numDOFs_1D)
1868 psi.resize(numDOFs_1D);
1869 for (
int i=0; i<NNZ_1D; i++)
1871 uStar_dMatrix[i]=0.;
1872 vStar_dMatrix[i]=0.;
1876 for (
int i=0; i<numDOFs_1D; i++)
1880 entropyResidualPerNode[i]=0.;
1890 double mesh_volume_conservation=0.0,
1891 mesh_volume_conservation_weak=0.0,
1892 mesh_volume_conservation_err_max=0.0,
1893 mesh_volume_conservation_err_max_weak=0.0;
1894 double globalConservationError=0.0;
1895 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
1896 for(
int eN=0;eN<nElements_global;eN++)
1898 double elementTransport[nDOF_test_element][nDOF_trial_element];
1899 double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1901 double elementResidual_p[nDOF_test_element],elementResidual_mesh[nDOF_test_element],
1902 elementResidual_u[nDOF_test_element],
1903 elementResidual_v[nDOF_test_element],
1904 mom_u_source_i[nDOF_test_element],
1905 mom_v_source_i[nDOF_test_element],
1906 betaDrag_i[nDOF_test_element],
1907 vos_i[nDOF_test_element],
1908 phisErrorElement[nDOF_test_element],
1910 elementEntropyResidual[nDOF_test_element],
1913 double element_active=1.0;
1914 double mesh_volume_conservation_element=0.0,
1915 mesh_volume_conservation_element_weak=0.0;
1917 double linVisc_eN = 0, nlinVisc_eN_num = 0, nlinVisc_eN_den = 0;
1919 double det_hess_uStar_Ke=0.0, det_hess_vStar_Ke=0.0, area_Ke=0.0;
1920 for (
int i=0;i<nDOF_test_element;i++)
1922 int eN_i = eN*nDOF_test_element+i;
1923 elementResidual_p_save[eN_i]=0.0;
1924 elementResidual_mesh[i]=0.0;
1925 elementResidual_p[i]=0.0;
1926 elementResidual_u[i]=0.0;
1927 elementResidual_v[i]=0.0;
1928 mom_u_source_i[i]=0.0;
1929 mom_v_source_i[i]=0.0;
1932 phisErrorElement[i]=0.0;
1934 elementEntropyResidual[i]=0.0;
1935 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1937 for (
int j=0;j<nDOF_trial_element;j++)
1939 elementTransport[i][j]=0.0;
1940 elementTransposeTransport[i][j]=0.0;
1945 if(use_ball_as_particle==1)
1947 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1949 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1950 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1951 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1952 phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]]);
1959 double _distance[nDOF_mesh_trial_element]={0.0};
1961 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1963 if(use_ball_as_particle==1)
1966 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1967 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1968 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1973 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
1975 if ( _distance[I] > 0)
1978 if (pos_counter == 3)
1980 element_active = 1.0;
1982 else if (pos_counter == 0)
1984 element_active = 1.0;
1988 element_active = 1.0;
1991 double sub_mesh_dof[6*3], sub_u_dof[15], sub_v_dof[15], sub_phi_dof[6], sub_p_dof[6];
1992 int boundaryNodes[6] = {0,0,0,0,0,0};
1993 std::vector<int> ls_nodes;
1994 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1996 for (
int K=0;K<nDOF_mesh_trial_element;K++)
2004 const double eps = 1.0e-4;
2005 double delta_phi=0.0,theta;
2006 delta_phi = _distance[(I+1)%3] - _distance[I];
2007 if (fabs(delta_phi) > eps)
2010 theta = -_distance[I]/delta_phi;
2011 if (theta > 1.0-eps || theta < eps)
2013 if (theta > 1.0-eps && theta <= 1.0)
2015 ls_nodes.push_back((I+1)%3);
2019 else if (theta > 0.0 && theta < eps)
2021 ls_nodes.push_back(I);
2029 boundaryNodes[3+I]=1;
2030 ls_nodes.push_back(3+I);
2036 if (fabs(_distance[I]) <= eps)
2039 boundaryNodes[3+I]=1;
2040 boundaryNodes[(I+1)%3]=1;
2041 ls_nodes.push_back(I);
2042 ls_nodes.push_back((I+1)%3);
2045 assert(theta <= 1.0);
2046 GI[3*3 + I*3 + I] = 1.0-theta;
2047 GI[3*3 + I*3 + (I+1)%3] = theta;
2048 GI[3*3 + I*3 + (I+2)%3] = 0.0;
2050 if (ls_nodes.size() != 2)
2052 std::cout<<
"level set nodes not 2 "<<ls_nodes.size()<<std::endl;
2053 for(
int i=0;i<ls_nodes.size();i++)
2054 std::cout<<ls_nodes[i]<<std::endl;
2055 std::sort(ls_nodes.begin(),ls_nodes.end());
2057 int sub_mesh_l2g[12] = {0,3,5,
2061 for (
int I=0; I<6; I++)
2063 sub_phi_dof[I] = 0.0;
2065 for (
int K=0; K<3; K++)
2066 sub_mesh_dof[I*3+K] = 0.0;
2067 for (
int J=0; J<3; J++)
2069 for (
int K=0; K<3; K++)
2071 sub_mesh_dof[I*3+K] += GI[I*3+J]*mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+J]+K];
2073 sub_phi_dof[I] += GI[I*3+J]*phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+J]];
2074 sub_p_dof[I] += GI[I*3+J]*p_dof[p_l2g[eN*nDOF_per_element_pressure+J]];
2077 int L = ls_nodes[0],
R=ls_nodes[1];
2078 double DX=sub_mesh_dof[
L*3+0] - sub_mesh_dof[
R*3+0];
2079 double DY=sub_mesh_dof[
L*3+1] - sub_mesh_dof[
R*3+1];
2080 double DS = std::sqrt(DX*DX+DY*DY);
2081 double nx = -DY/DS, ny = DX/DS;
2082 double nxL,nyL,nxR,nyR;
2085 sub_mesh_dof[
L*3+0],sub_mesh_dof[
L*3+1],0.0,
2089 sub_mesh_dof[
R*3+0],sub_mesh_dof[
R*3+1],0.0,
2093 double n_fluid_sign = -
sgn(nx*0.5*(nxL+nxR)+ny*0.5*(nyL+nyR));
2098 cut_cell_boundary_length += DS;
2099 p_force_x += sub_p_dof[
L]*nx*0.5*DS + sub_p_dof[
R]*nx*0.5*DS;
2100 p_force_y += sub_p_dof[
L]*ny*0.5*DS + sub_p_dof[
R]*ny*0.5*DS;
2107 double lagrangeNodes[9*3];
2108 for (
int K=0;K<3;K++)
2110 lagrangeNodes[0*3+K] = 0.5*(sub_mesh_dof[0*3+K] + sub_mesh_dof[3*3+0*3+K]);
2111 lagrangeNodes[1*3+K] = 0.5*(sub_mesh_dof[1*3+K] + sub_mesh_dof[3*3+0*3+K]);
2112 lagrangeNodes[2*3+K] = 0.5*(sub_mesh_dof[1*3+K] + sub_mesh_dof[3*3+1*3+K]);
2113 lagrangeNodes[3*3+K] = 0.5*(sub_mesh_dof[2*3+K] + sub_mesh_dof[3*3+1*3+K]);
2114 lagrangeNodes[4*3+K] = 0.5*(sub_mesh_dof[2*3+K] + sub_mesh_dof[3*3+2*3+K]);
2115 lagrangeNodes[5*3+K] = 0.5*(sub_mesh_dof[0*3+K] + sub_mesh_dof[3*3+2*3+K]);
2116 lagrangeNodes[6*3+K] = 0.5*(sub_mesh_dof[3*3+0*3+K] + sub_mesh_dof[3*3+1*3+K]);
2117 lagrangeNodes[7*3+K] = 0.5*(sub_mesh_dof[3*3+1*3+K] + sub_mesh_dof[3*3+2*3+K]);
2118 lagrangeNodes[8*3+K] = 0.5*(sub_mesh_dof[3*3+2*3+K] + sub_mesh_dof[3*3+0*3+K]);
2121 for (
int I=0;I<6;I++)
2123 baryCoords(&sub_mesh_dof[0],&sub_mesh_dof[1*3],&sub_mesh_dof[2*3],&sub_mesh_dof[I*3],lambda);
2125 G2I[I*6+0] = lambda[0]*(2.0*lambda[0] - 1.0);
2126 G2I[I*6+1] = lambda[1]*(2.0*lambda[1] - 1.0);
2127 G2I[I*6+2] = lambda[2]*(2.0*lambda[2] - 1.0);
2128 G2I[I*6+3] = 4.0*lambda[0]*lambda[1];
2129 G2I[I*6+4] = 4.0*lambda[1]*lambda[2];
2130 G2I[I*6+5] = 4.0*lambda[2]*lambda[0];
2132 for (
int I=0;I<9;I++)
2134 baryCoords(&sub_mesh_dof[0],&sub_mesh_dof[1*3],&sub_mesh_dof[2*3],&lagrangeNodes[I*3],lambda);
2135 G2I[6*6 + I*6 + 0] = lambda[0]*(2.0*lambda[0] - 1.0);
2136 G2I[6*6 + I*6 + 1] = lambda[1]*(2.0*lambda[1] - 1.0);
2137 G2I[6*6 + I*6 + 2] = lambda[2]*(2.0*lambda[2] - 1.0);
2138 G2I[6*6 + I*6 + 3] = 4.0*lambda[0]*lambda[1];
2139 G2I[6*6 + I*6 + 4] = 4.0*lambda[1]*lambda[2];
2140 G2I[6*6 + I*6 + 5] = 4.0*lambda[2]*lambda[0];
2142 for (
int I=0; I<15; I++)
2146 for (
int J=0; J<6; J++)
2148 sub_u_dof[I] += G2I[I*6+J]*u_dof[vel_l2g[eN*nDOF_trial_element+J]];
2149 sub_v_dof[I] += G2I[I*6+J]*v_dof[vel_l2g[eN*nDOF_trial_element+J]];
2152 for (
int esN=0;esN<4;esN++)
2154 std::cout<<sub_mesh_l2g[esN*3]<<
'\t'<<sub_mesh_l2g[esN*3+1]<<
'\t'<<sub_mesh_l2g[esN*3+2]<<std::endl;
2156 for (
int I=0; I<6; I++)
2158 std::cout<<sub_mesh_dof[I*3+0]<<
'\t'<<sub_mesh_dof[I*3+1]<<
'\t'<<sub_mesh_dof[I*3+2]<<
'\t'<<boundaryNodes[I]<<
'\t'<<sub_phi_dof[I]<<
'\t'<<sub_p_dof[I]<<
'\t'<<sub_u_dof[I]<<
'\t'<<sub_v_dof[I]<<
'\t'<<G2I[I*6+0]<<
'\t'<<G2I[I*6+1]<<
'\t'<<G2I[I*6+2]<<
'\t'<<G2I[I*6+3]<<
'\t'<<G2I[I*6+4]<<
'\t'<<G2I[I*6+5]<<std::endl;
2167 double _distance[nDOF_mesh_trial_element]={0.0};
2169 for (
int I=0;I<nDOF_mesh_trial_element;I++)
2171 if(use_ball_as_particle==1)
2174 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
2175 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
2176 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
2181 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
2183 if ( _distance[I] >= 0)
2186 if (pos_counter == 2)
2190 for (
int I=0;I<nDOF_mesh_trial_element;I++)
2193 if (_distance[I] < 0)
2199 assert(opp_node >=0);
2200 assert(opp_node <nDOF_mesh_trial_element);
2205 const int ebN = elementBoundariesArray[eN*nDOF_mesh_trial_element+opp_node];
2206 const int eN_oppo = (eN == elementBoundaryElementsArray[ebN*2+0])?elementBoundaryElementsArray[ebN*2+1]:elementBoundaryElementsArray[ebN*2+0];
2207 if((mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%3]<nNodes_owned
2208 || mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%3]<nNodes_owned)
2214 if (eN == elementBoundaryElementsArray[ebN*2+0])
2221 if(use_ball_as_particle==1)
2223 double middle_point_coord[3]={0.0};
2224 double middle_point_distance;
2225 middle_point_coord[0] = 0.5*(mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%3]+0]+mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%3]+0]);
2226 middle_point_coord[1] = 0.5*(mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%3]+1]+mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%3]+1]);
2228 middle_point_coord[0],middle_point_coord[1],middle_point_coord[2],
2229 middle_point_distance);
2239 double distance=1e10, distance_to_ith_particle;
2240 for (
int i=0;i<nParticles;++i)
2242 distance_to_ith_particle=particle_signed_distances[i*nElements_global*nQuadraturePoints_element
2243 +eN*nQuadraturePoints_element
2245 if (distance_to_ith_particle<distance)
2247 distance = distance_to_ith_particle;
2257 if(ebN<nElementBoundaries_owned)
2259 assert(eN_oppo==-1);
2263 else if (pos_counter == 3)
2266 for (
int i=0;i<nDOF_test_element;i++)
2268 isActiveDOF[offset_u+stride_u*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
2269 isActiveDOF[offset_v+stride_v*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
2277 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
2278 for (
int j=0;j<nDOF_mesh_trial_element;j++)
2280 int eN_j = eN*nDOF_mesh_trial_element+j;
2281 element_phi[j] = phi_dof[p_l2g[eN_j]];
2282 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
2284 double element_nodes[nDOF_mesh_trial_element*3];
2285 for (
int i=0;i<nDOF_mesh_trial_element;i++)
2287 int eN_i=eN*nDOF_mesh_trial_element+i;
2288 for(
int I=0;I<3;I++)
2289 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
2291 gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
2292 gf.
calculate(element_phi, element_nodes, x_ref.data(),
false);
2296 for(
int k=0;k<nQuadraturePoints_element;k++)
2301 int eN_k = eN*nQuadraturePoints_element+k,
2302 eN_k_nSpace = eN_k*nSpace,
2304 eN_nDOF_trial_element = eN*nDOF_trial_element;
2305 double p=0.0,
u=0.0,
v=0.0,
w=0.0,un=0.0,vn=0.0,wn=0.0,
2306 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
2315 dmass_adv_u[nSpace],
2316 dmass_adv_v[nSpace],
2317 dmass_adv_w[nSpace],
2319 dmom_u_adv_u[nSpace],
2320 dmom_u_adv_v[nSpace],
2321 dmom_u_adv_w[nSpace],
2323 dmom_v_adv_u[nSpace],
2324 dmom_v_adv_v[nSpace],
2325 dmom_v_adv_w[nSpace],
2327 dmom_w_adv_u[nSpace],
2328 dmom_w_adv_v[nSpace],
2329 dmom_w_adv_w[nSpace],
2330 mom_uu_diff_ten[nSpace],
2331 mom_vv_diff_ten[nSpace],
2332 mom_ww_diff_ten[nSpace],
2343 dmom_u_ham_grad_p[nSpace],
2344 dmom_u_ham_grad_u[nSpace],
2346 dmom_v_ham_grad_p[nSpace],
2347 dmom_v_ham_grad_v[nSpace],
2349 dmom_w_ham_grad_p[nSpace],
2350 dmom_w_ham_grad_w[nSpace],
2361 Lstar_u_p[nDOF_test_element],
2362 Lstar_v_p[nDOF_test_element],
2363 Lstar_w_p[nDOF_test_element],
2364 Lstar_u_u[nDOF_test_element],
2365 Lstar_v_v[nDOF_test_element],
2366 Lstar_w_w[nDOF_test_element],
2367 Lstar_p_u[nDOF_test_element],
2368 Lstar_p_v[nDOF_test_element],
2369 Lstar_p_w[nDOF_test_element],
2374 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
2375 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2378 jacInv[nSpace*nSpace],
2379 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
2380 vel_hess_trial[nDOF_trial_element*
nSpace2],
2381 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
2382 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
2383 u_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2384 v_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2390 dmom_u_source[nSpace],
2391 dmom_v_source[nSpace],
2392 dmom_w_source[nSpace],
2396 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
2398 ck.calculateMapping_element(eN,
2402 mesh_trial_ref.data(),
2403 mesh_grad_trial_ref.data(),
2408 ck.calculateH_element(eN,
2410 nodeDiametersArray.data(),
2412 mesh_trial_ref.data(),
2414 ck.calculateMappingVelocity_element(eN,
2416 mesh_velocity_dof.data(),
2418 mesh_trial_ref.data(),
2423 dV = fabs(jacDet)*dV_ref[k];
2424 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2427 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2428 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2429 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2433 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
2434 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
2439 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
2440 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
2443 ck.valFromDOF(u_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],un);
2444 ck.valFromDOF(v_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],vn);
2448 for (
int I=0;I<nSpace;I++)
2449 grad_p[I] = q_grad_p[eN_k_nSpace + I];
2450 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
2451 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
2452 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
2453 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
2454 ck.hessFromDOF(uStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_uStar);
2455 ck.hessFromDOF(vStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_vStar);
2458 for (
int j=0;j<nDOF_trial_element;j++)
2461 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
2462 for (
int I=0;I<nSpace;I++)
2465 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
2466 if (ARTIFICIAL_VISCOSITY==4)
2469 u_times_vel_grad_test_dV[j*nSpace+I] =
2470 u*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_u[I];
2471 v_times_vel_grad_test_dV[j*nSpace+I] =
2472 v*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_v[I];
2479 if (ARTIFICIAL_VISCOSITY==3)
2481 det_hess_uStar_Ke += (hess_uStar[0]*hess_uStar[3] - hess_uStar[2]*hess_uStar[1])*dV;
2482 det_hess_vStar_Ke += (hess_vStar[0]*hess_vStar[3] - hess_vStar[2]*hess_vStar[1])*dV;
2486 double div_mesh_velocity=0.0;
2487 int NDOF_MESH_TRIAL_ELEMENT=3;
2488 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
2490 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
2491 div_mesh_velocity +=
2492 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*nSpace+0] +
2493 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*nSpace+1];
2495 mesh_volume_conservation_element += (alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity)*dV;
2496 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
2498 porosity = 1.0 - q_vos[eN_k];
2507 double distance_to_omega_solid = 1e10;
2508 if(use_ball_as_particle==1)
2512 distance_to_omega_solid);
2516 for (
int i = 0; i < nParticles; i++)
2518 double distance_to_i_th_solid = particle_signed_distances[i * nElements_global * nQuadraturePoints_element + eN_k];
2519 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
2522 phi_solid[eN_k] = distance_to_omega_solid;
2534 elementDiameter[eN],
2535 smagorinskyConstant,
2536 turbulenceClosureModel,
2541 &normal_phi[eN_k_nSpace],
2542 distance_to_omega_solid,
2555 q_velocity_sge[eN_k_nSpace+0],
2556 q_velocity_sge[eN_k_nSpace+1],
2557 q_velocity_sge[eN_k_nSpace+1],
2558 q_eddy_viscosity[eN_k],
2605 MULTIPLY_EXTERNAL_FORCE_BY_DENSITY,
2609 MATERIAL_PARAMETERS_AS_FUNCTION,
2610 density_as_function[eN_k],
2611 dynamic_viscosity_as_function[eN_k],
2614 use_ball_as_particle,
2617 ball_velocity.data(),
2618 ball_angular_velocity.data(),
2619 INT_BY_PARTS_PRESSURE);
2622 mass_source = q_mass_source[eN_k];
2623 for (
int I=0;I<nSpace;I++)
2625 dmom_u_source[I] = 0.0;
2626 dmom_v_source[I] = 0.0;
2627 dmom_w_source[I] = 0.0;
2638 q_eddy_viscosity[eN_k],
2645 q_velocity_sge[eN_k_nSpace+0],
2646 q_velocity_sge[eN_k_nSpace+1],
2647 q_velocity_sge[eN_k_nSpace+1],
2648 eps_solid[elementFlags[eN]],
2650 q_velocity_solid[eN_k_nSpace+0],
2651 q_velocity_solid[eN_k_nSpace+1],
2652 q_velocity_solid[eN_k_nSpace+1],
2653 q_velocityStar_solid[eN_k_nSpace+0],
2654 q_velocityStar_solid[eN_k_nSpace+1],
2655 q_velocityStar_solid[eN_k_nSpace+1],
2662 q_grad_vos[eN_k_nSpace+0],
2663 q_grad_vos[eN_k_nSpace+1],
2664 q_grad_vos[eN_k_nSpace+1]);
2665 double C_particles=0.0;
2666 if(nParticles > 0 && USE_SBM==0)
2671 nQuadraturePoints_global,
2672 &particle_signed_distances[eN_k],
2673 &particle_signed_distance_normals[eN_k_3d],
2674 &particle_velocities[eN_k_3d],
2675 particle_centroids.data(),
2676 use_ball_as_particle,
2679 ball_velocity.data(),
2680 ball_angular_velocity.data(),
2682 particle_penalty_constant/h_phi,
2683 particle_alpha/h_phi,
2684 particle_beta/h_phi,
2701 q_velocity_sge[eN_k_nSpace+0],
2702 q_velocity_sge[eN_k_nSpace+1],
2703 q_velocity_sge[eN_k_nSpace+1],
2726 particle_netForces.data(),
2727 particle_netMoments.data(),
2728 particle_surfaceArea.data());
2733 nQuadraturePoints_global,
2734 &particle_signed_distances[eN_k],
2735 &particle_signed_distance_normals[eN_k_3d],
2736 &particle_velocities[eN_k_3d],
2737 particle_centroids.data(),
2738 use_ball_as_particle,
2741 ball_velocity.data(),
2742 ball_angular_velocity.data(),
2743 particle_penalty_constant/h_phi,
2744 particle_alpha/h_phi,
2745 particle_beta/h_phi,
2762 q_velocity_sge[eN_k_nSpace+0],
2763 q_velocity_sge[eN_k_nSpace+1],
2764 q_velocity_sge[eN_k_nSpace+1],
2769 particle_netForces.data(),
2770 particle_netMoments.data());
2772 if (turbulenceClosureModel >= 3)
2774 const double c_mu = 0.09;
2789 &q_turb_var_grad_0[eN_k_nSpace],
2790 q_eddy_viscosity[eN_k],
2808 q_mom_u_acc[eN_k] = mom_u_acc;
2809 q_mom_v_acc[eN_k] = mom_v_acc;
2812 q_mass_adv[eN_k_nSpace+0] =
u;
2813 q_mass_adv[eN_k_nSpace+1] =
v;
2818 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
2819 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
2821 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
2822 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
2825 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
2826 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
2828 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
2829 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
2841 if (q_dV_last[eN_k] <= -100)
2842 q_dV_last[eN_k] = dV;
2845 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2851 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2865 mom_u_acc_t *= dmom_u_acc_u;
2866 mom_v_acc_t *= dmom_v_acc_v;
2872 ck.Mass_strong(-q_dvos_dt[eN_k]) +
2873 ck.Advection_strong(dmass_adv_u,grad_u) +
2874 ck.Advection_strong(dmass_adv_v,grad_v) +
2876 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
2878 ck.Reaction_strong(mass_source);
2881 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
2882 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
2886 ck.Mass_strong(mom_u_acc_t) +
2887 ck.Advection_strong(dmom_adv_sge,grad_u) +
2888 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
2889 ck.Reaction_strong(mom_u_source) -
2890 ck.Reaction_strong(
u*div_mesh_velocity);
2893 ck.Mass_strong(mom_v_acc_t) +
2894 ck.Advection_strong(dmom_adv_sge,grad_v) +
2895 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
2896 ck.Reaction_strong(mom_v_source) -
2897 ck.Reaction_strong(
v*div_mesh_velocity);
2907 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
2909 elementDiameter[eN],
2914 dmom_u_ham_grad_p[0],
2924 dmom_u_ham_grad_p[0],
2929 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
2930 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
2943 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
2944 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
2947 mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge[eN_k_nSpace+0]);
2948 mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge[eN_k_nSpace+0]);
2952 for (
int i=0;i<nDOF_test_element;i++)
2954 int i_nSpace = i*nSpace;
2959 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2960 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2962 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2963 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2967 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2968 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2973 if (ARTIFICIAL_VISCOSITY==0 || ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2975 q_numDiff_u[eN_k] = 0;
2976 q_numDiff_v[eN_k] = 0;
2977 q_numDiff_w[eN_k] = 0;
2979 else if (ARTIFICIAL_VISCOSITY==1)
2981 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v);
2982 q_numDiff_u[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2983 (1.0-useMetrics)*hFactor*hFactor*elementDiameter[eN]*elementDiameter[eN]);
2984 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
2985 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
2989 double rho = q_rho[eN_k];
2990 double mu = q_rho[eN_k]*q_nu[eN_k];
2992 double vel2 =
u*
u +
v*
v;
2996 porosity*rho*((
u-un)/dt + (
u*grad_u[0]+
v*grad_u[1]) - g[0])
2997 + (KILL_PRESSURE_TERM == 1 ? 0. : 1.)*grad_p[0]
2998 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcex[eN_k]
2999 - mu*(hess_u[0] + hess_u[3])
3000 - mu*(hess_u[0] + hess_v[2]);
3003 porosity*rho*((
v-vn)/dt + (
u*grad_v[0]+
v*grad_v[1]) - g[1])
3004 + (KILL_PRESSURE_TERM == 1 ? 0. : 1.)*grad_p[1]
3005 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcey[eN_k]
3006 - mu*(hess_v[0] + hess_v[3])
3007 - mu*(hess_u[1] + hess_v[3]);
3010 double entRes_times_u = Res_in_x*
u + Res_in_y*
v;
3012 double hK = elementDiameter[eN]/order_polynomial;
3013 q_numDiff_u[eN_k] = fmin(
cMax*porosity*rho*hK*std::sqrt(vel2),
3014 cE*hK*hK*fabs(entRes_times_u)/(vel2+1E-10));
3015 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
3016 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
3020 linVisc_eN = fmax(porosity*rho*std::sqrt(vel2),linVisc_eN);
3021 nlinVisc_eN_num = fmax(fabs(entRes_times_u),nlinVisc_eN_num);
3022 nlinVisc_eN_den = fmax(vel2,nlinVisc_eN_den);
3033 q_velocity[eN_k_nSpace+0]=
u;
3034 q_velocity[eN_k_nSpace+1]=
v;
3036 for (
int I=0;I<nSpace;I++)
3038 q_grad_u[eN_k_nSpace+I] = grad_u[I];
3039 q_grad_v[eN_k_nSpace+I] = grad_v[I];
3043 q_divU[eN_k] = q_grad_u[eN_k_nSpace+0] + q_grad_v[eN_k_nSpace+1];
3046 double unit_normal[nSpace];
3047 double norm_grad_phi = 0.;
3048 for (
int I=0;I<nSpace;I++)
3049 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
3050 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
3051 for (
int I=0;I<nSpace;I++)
3052 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
3056 v1[0]=1.-unit_normal[0]*unit_normal[0];
3057 v1[1]=-unit_normal[0]*unit_normal[1];
3060 v2[0]=-unit_normal[0]*unit_normal[1];
3061 v2[1]=1.-unit_normal[1]*unit_normal[1];
3062 double delta =
gf.
D(eps_mu,
phi[eN_k]);
3063 double vel_tgrad_test_i[nSpace], tgrad_u[nSpace], tgrad_v[nSpace];
3072 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3074 velStar[0] = q_velocity_sge[eN_k_nSpace+0];
3075 velStar[1] = q_velocity_sge[eN_k_nSpace+1];
3078 for(
int i=0;i<nDOF_test_element;i++)
3080 int i_nSpace=i*nSpace;
3082 &vel_grad_trial[i_nSpace],
3084 phisErrorElement[i]+=std::abs(phisError[eN_k_nSpace+0])*p_test_dV[i];
3102 elementResidual_u[i] +=
3103 ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
3104 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
3105 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]) +
3106 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]) +
3108 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
3109 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
3110 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+0] : 0.) +
3112 USE_SUPG*
ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
3113 ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]) +
3115 ck.NumericalDiffusion(delta*sigma*dV,v1,vel_tgrad_test_i) +
3116 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_u,vel_tgrad_test_i);
3117 mom_u_source_i[i] +=
ck.Reaction_weak(mom_u_source,vel_test_dV[i]);
3118 betaDrag_i[i] +=
ck.Reaction_weak(dmom_u_source[0],
3120 vos_i[i] +=
ck.Reaction_weak(1.0-porosity,
3123 elementResidual_v[i] +=
3124 ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
3125 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
3126 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]) +
3127 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]) +
3129 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
3130 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
3131 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+1] : 0.) +
3133 USE_SUPG*
ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
3134 ck.NumericalDiffusion(q_numDiff_v_last[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]) +
3136 ck.NumericalDiffusion(delta*sigma*dV,v2,vel_tgrad_test_i) +
3137 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_v,vel_tgrad_test_i);
3138 mom_v_source_i[i] +=
ck.Reaction_weak(mom_v_source,vel_test_dV[i]);
3152 if (ARTIFICIAL_VISCOSITY==4)
3156 elementEntropyResidual[i] +=
3158 ck.Mass_weak(mom_u_acc_t,
u*vel_test_dV[i]) +
3159 ck.Advection_weak(mom_u_adv,&u_times_vel_grad_test_dV[i_nSpace])+
3160 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),
3161 sdInfo_u_u_colind.data(),
3164 &u_times_vel_grad_test_dV[i_nSpace]) +
3165 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),
3166 sdInfo_u_v_colind.data(),
3169 &u_times_vel_grad_test_dV[i_nSpace]) +
3170 ck.Reaction_weak(mom_u_source,
u*vel_test_dV[i]) +
3171 ck.Hamiltonian_weak(mom_u_ham,
u*vel_test_dV[i])
3173 ck.Mass_weak(mom_v_acc_t,
v*vel_test_dV[i]) +
3174 ck.Advection_weak(mom_v_adv,&v_times_vel_grad_test_dV[i_nSpace])+
3175 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),
3176 sdInfo_v_u_colind.data(),
3179 &v_times_vel_grad_test_dV[i_nSpace])+
3180 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),
3181 sdInfo_v_v_colind.data(),
3184 &v_times_vel_grad_test_dV[i_nSpace])+
3185 ck.Reaction_weak(mom_v_source,
v*vel_test_dV[i]) +
3186 ck.Hamiltonian_weak(mom_v_ham,
v*vel_test_dV[i]);
3188 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3190 for(
int j=0;j<nDOF_trial_element;j++)
3192 int j_nSpace = j*nSpace;
3193 int i_nSpace = i*nSpace;
3194 elementTransport[i][j] +=
3195 q_rho[eN_k]*porosity*
3196 ck.AdvectionJacobian_strong(velStar,
3197 &vel_grad_test_dV[j_nSpace])
3198 *vel_trial_ref[k*nDOF_trial_element+i];
3199 elementTransposeTransport[i][j] +=
3200 q_rho[eN_k]*porosity*
3201 ck.AdvectionJacobian_strong(velStar,
3202 &vel_grad_test_dV[i_nSpace])
3203 *vel_trial_ref[k*nDOF_trial_element+j];
3208 element_uStar_He[eN] = det_hess_uStar_Ke/area_Ke;
3209 element_vStar_He[eN] = det_hess_vStar_Ke/area_Ke;
3214 double hK = elementDiameter[eN];
3215 double artVisc = fmin(
cMax*hK*linVisc_eN,
3216 cE*hK*hK*nlinVisc_eN_num/(nlinVisc_eN_den+1E-10));
3217 for(
int k=0;k<nQuadraturePoints_element;k++)
3219 int eN_k = eN*nQuadraturePoints_element+k;
3220 q_numDiff_u[eN_k] = artVisc;
3221 q_numDiff_v[eN_k] = artVisc;
3222 q_numDiff_w[eN_k] = artVisc;
3228 for(
int i=0;i<nDOF_test_element;i++)
3230 int eN_i=eN*nDOF_test_element+i;
3231 phisErrorNodal[vel_l2g[eN_i]]+= element_active*phisErrorElement[i];
3235 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=element_active*elementResidual_u[i];
3236 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=element_active*elementResidual_v[i];
3238 ncDrag[offset_u+stride_u*vel_l2g[eN_i]]+=mom_u_source_i[i];
3239 ncDrag[offset_v+stride_v*vel_l2g[eN_i]]+=mom_v_source_i[i];
3240 betaDrag[vel_l2g[eN_i]] += betaDrag_i[i];
3241 vos_vel_nodes[vel_l2g[eN_i]] += vos_i[i];
3244 if (ARTIFICIAL_VISCOSITY==3)
3246 uStar_hi[vel_l2g[eN_i]] += element_uStar_He[eN];
3247 vStar_hi[vel_l2g[eN_i]] += element_vStar_He[eN];
3248 den_hi[vel_l2g[eN_i]] += 1;
3250 if (ARTIFICIAL_VISCOSITY==4)
3253 entropyResidualPerNode[vel_l2g[eN_i]] += elementEntropyResidual[i];
3255 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3257 for (
int j=0;j<nDOF_trial_element;j++)
3259 int eN_i_j = eN_i*nDOF_trial_element+j;
3261 + csrColumnOffsets_1D[eN_i_j]]
3262 += elementTransport[i][j];
3265 + csrColumnOffsets_1D[eN_i_j]]
3266 += elementTransposeTransport[i][j];
3277 std::cout<<std::flush;
3279 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3282 for (
int i=0; i<numDOFs_1D; i++)
3284 if (ARTIFICIAL_VISCOSITY==4)
3287 double max_u2i = (std::pow(u_dof[i],2.) +
3288 std::pow(v_dof[i],2.));
3289 double min_u2i = max_u2i;
3290 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3292 int j = colind_1D[offset];
3293 double u2j = (std::pow(u_dof[j],2.) +
3294 std::pow(v_dof[j],2.));
3295 max_u2i = fmax(max_u2i,u2j);
3296 min_u2i = fmin(min_u2i,u2j);
3298 double normi = 0.5*(max_u2i + min_u2i) + 1E-10;
3299 entropyResidualPerNode[i] = fabs(entropyResidualPerNode[i])/normi;
3304 double uStari = uStar_dof[i];
3305 double vStari = vStar_dof[i];
3307 double u_beta_numerator = 0., u_beta_denominator = 0.;
3308 double v_beta_numerator = 0., v_beta_denominator = 0.;
3311 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3313 int j = colind_1D[offset];
3314 double uStarj = uStar_dof[j];
3315 double vStarj = vStar_dof[j];
3318 u_beta_numerator += (uStarj - uStari);
3319 u_beta_denominator += fabs(uStarj - uStari);
3321 v_beta_numerator += (vStarj - vStari);
3322 v_beta_denominator += fabs(vStarj - vStari);
3324 double u_beta = fabs(u_beta_numerator)/(u_beta_denominator+1E-10);
3325 double v_beta = fabs(v_beta_numerator)/(v_beta_denominator+1E-10);
3344 if (ARTIFICIAL_VISCOSITY==3)
3346 for(
int eN=0;eN<nElements_global;eN++)
3348 double uStar_He = element_uStar_He[eN];
3349 double vStar_He = element_vStar_He[eN];
3350 for(
int i=0;i<nDOF_test_element;i++)
3352 int eN_i=eN*nDOF_test_element+i;
3353 int gi = vel_l2g[eN_i];
3361 if (ARTIFICIAL_VISCOSITY==3)
3363 for (
int i=0; i<numDOFs_1D; i++)
3368 if (isBoundary_1D[i] == 1)
3392 for (
int i=0; i<numDOFs_1D; i++)
3395 double uStar_dii = 0;
3396 double vStar_dii = 0;
3397 double ui = u_dof[i];
3398 double vi = v_dof[i];
3400 double ith_u_dissipative_term = 0;
3401 double ith_v_dissipative_term = 0;
3406 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3408 int j = colind_1D[offset];
3411 double uj = u_dof[j];
3412 double vj = v_dof[j];
3417 if (ARTIFICIAL_VISCOSITY==4)
3419 double dEVij = fmax(laggedEntropyResidualPerNode[i],
3420 laggedEntropyResidualPerNode[j]);
3423 uStar_dMatrix[ij] = fmin(dLij,
cE*dEVij);
3424 vStar_dMatrix[i] = uStar_dMatrix[ij];
3433 uStar_dii -= uStar_dMatrix[ij];
3434 vStar_dii -= vStar_dMatrix[ij];
3436 ith_u_dissipative_term += uStar_dMatrix[ij]*(uj-ui);
3437 ith_v_dissipative_term += vStar_dMatrix[ij]*(vj-vi);
3446 uStar_dMatrix[ii] = uStar_dii;
3447 vStar_dMatrix[ii] = vStar_dii;
3448 globalResidual[offset_u+stride_u*i] += -ith_u_dissipative_term;
3449 globalResidual[offset_v+stride_v*i] += -ith_v_dissipative_term;
3460 std::memset(particle_netForces.data(),0,nParticles*3*
sizeof(
double));
3461 std::memset(particle_netMoments.data(),0,nParticles*3*
sizeof(
double));
3466 double Fx = 0.0, Fy = 0.0, Fxp = 0.0, Fyp = 0.0, surfaceArea=0.0, Mz = 0.0;
3470 eN_nDOF_trial_element = eN*nDOF_trial_element;
3471 double elementResidual_mesh[nDOF_test_element],
3472 elementResidual_p[nDOF_test_element],
3473 elementResidual_u[nDOF_test_element],
3474 elementResidual_v[nDOF_test_element],
3481 for (
int i=0;i<nDOF_test_element;i++)
3483 elementResidual_mesh[i]=0.0;
3484 elementResidual_p[i]=0.0;
3485 elementResidual_u[i]=0.0;
3486 elementResidual_v[i]=0.0;
3489 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3491 int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
3493 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3494 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3502 jac_ext[nSpace*nSpace],
3504 jacInv_ext[nSpace*nSpace],
3505 boundaryJac[nSpace*(nSpace-1)],
3506 metricTensor[(nSpace-1)*(nSpace-1)],
3507 metricTensorDetSqrt,
3508 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
3509 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
3510 vel_grad_test_dS[nDOF_trial_element*nSpace],
3511 normal[2],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3512 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3513 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;
3515 ck.calculateMapping_elementBoundary(eN,
3521 mesh_trial_trace_ref.data(),
3522 mesh_grad_trial_trace_ref.data(),
3523 boundaryJac_ref.data(),
3529 metricTensorDetSqrt,
3533 ck.calculateMappingVelocity_elementBoundary(eN,
3537 mesh_velocity_dof.data(),
3539 mesh_trial_trace_ref.data(),
3540 xt_ext,yt_ext,zt_ext,
3545 dS = metricTensorDetSqrt*dS_ref[kb];
3547 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3550 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3552 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3553 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3555 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3556 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3558 for (
int j=0;j<nDOF_trial_element;j++)
3560 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3561 for (
int I=0;I<nSpace;I++)
3562 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3566 double distance[2], P_normal[2], P_tangent[2];
3570 if(use_ball_as_particle==1)
3578 P_normal[0],P_normal[1]);
3580 ball_velocity.data(),ball_angular_velocity.data(),
3582 x_ext-dist*P_normal[0],
3583 y_ext-dist*P_normal[1],
3589 dist = ebq_global_phi_solid[ebN_kb];
3590 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*3+0];
3591 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*3+1];
3592 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*3+0];
3593 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*3+1];
3597 ck.calculateGScale(G,normal,h_penalty);
3601 assert(h_penalty>0.0);
3602 if (h_penalty < std::abs(dist))
3603 h_penalty = std::abs(dist);
3604 distance[0] = -P_normal[0]*dist;
3605 distance[1] = -P_normal[1]*dist;
3606 P_tangent[0] = -P_normal[1];
3607 P_tangent[1] = P_normal[0];
3609 double C_adim =
C_sbm*visco/h_penalty;
3610 double beta_adim =
beta_sbm*visco/h_penalty;
3615 const double u_m_uD[2] = {u_ext - bc_u_ext,v_ext - bc_v_ext};
3616 const double zero_vec[2]={0.,0.};
3619 for (
int i=0;i<nDOF_test_element;i++)
3621 int eN_i = eN*nDOF_test_element+i;
3623 int GlobPos_u = offset_u+stride_u*vel_l2g[eN_i];
3624 int GlobPos_v = offset_v+stride_v*vel_l2g[eN_i];
3625 double phi_i = vel_test_dS[i];
3626 double Gxphi_i = vel_grad_test_dS[i*nSpace+0];
3627 double Gyphi_i = vel_grad_test_dS[i*nSpace+1];
3628 double *grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
3630 const double grad_phi_i_dot_t =
get_dot_product(P_tangent,grad_phi_i);
3633 globalResidual[GlobPos_u] += C_adim*phi_i*u_m_uD[0];
3634 globalResidual[GlobPos_v] += C_adim*phi_i*u_m_uD[1];
3635 Fx += C_adim*phi_i*u_m_uD[0];
3636 Fy += C_adim*phi_i*u_m_uD[1];
3640 globalResidual[GlobPos_u] -= visco * phi_i*res[0];
3641 globalResidual[GlobPos_v] -= visco * phi_i*res[1];
3642 Fx -= visco * phi_i*res[0];
3643 Fy -= visco * phi_i*res[1];
3656 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*u_m_uD[0];
3657 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*u_m_uD[1];
3658 Fx += C_adim*grad_phi_i_dot_d*u_m_uD[0];
3659 Fy += C_adim*grad_phi_i_dot_d*u_m_uD[1];
3662 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*grad_u_d[0];
3663 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*grad_u_d[1];
3664 Fx += C_adim*grad_phi_i_dot_d*grad_u_d[0];
3665 Fy += C_adim*grad_phi_i_dot_d*grad_u_d[1];
3668 globalResidual[GlobPos_u] += C_adim*phi_i*grad_u_d[0];
3669 globalResidual[GlobPos_v] += C_adim*phi_i*grad_u_d[1];
3670 Fx += C_adim*phi_i*grad_u_d[0];
3671 Fy += C_adim*phi_i*grad_u_d[1];
3685 globalResidual[GlobPos_u] += beta_adim*grad_u_t[0]*grad_phi_i_dot_t;
3686 globalResidual[GlobPos_v] += beta_adim*grad_u_t[1]*grad_phi_i_dot_t;
3687 Fx += beta_adim*grad_u_t[0]*grad_phi_i_dot_t;
3688 Fy += beta_adim*grad_u_t[1]*grad_phi_i_dot_t;
3697 for (
int i=0; i<nDOF_per_element_pressure;++i)
3699 p_ext += p_dof[p_l2g[eN*nDOF_per_element_pressure+i]]*p_trial_trace_ref[ebN_local_kb*nDOF_per_element_pressure+i];
3701 double nx = P_normal[0];
3702 double ny = P_normal[1];
3708 if(use_ball_as_particle==1)
3718 Mz += r_x*Fy-r_y*Fx;
3721 && ebN < nElementBoundaries_owned)
3743 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3745 int ebN = exteriorElementBoundariesArray[ebNE],
3746 eN = elementBoundaryElementsArray[ebN*2+0],
3747 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
3748 eN_nDOF_trial_element = eN*nDOF_trial_element;
3749 double elementResidual_mesh[nDOF_test_element],
3750 elementResidual_p[nDOF_test_element],
3751 elementResidual_u[nDOF_test_element],
3752 elementResidual_v[nDOF_test_element],
3755 const double* elementResidual_w(NULL);
3756 for (
int i=0;i<nDOF_test_element;i++)
3758 elementResidual_mesh[i]=0.0;
3759 elementResidual_p[i]=0.0;
3760 elementResidual_u[i]=0.0;
3761 elementResidual_v[i]=0.0;
3764 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3766 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3767 ebNE_kb_nSpace = ebNE_kb*nSpace,
3768 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3769 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3779 dmom_u_acc_u_ext=0.0,
3781 dmom_v_acc_v_ext=0.0,
3783 dmom_w_acc_w_ext=0.0,
3784 mass_adv_ext[nSpace],
3785 dmass_adv_u_ext[nSpace],
3786 dmass_adv_v_ext[nSpace],
3787 dmass_adv_w_ext[nSpace],
3788 mom_u_adv_ext[nSpace],
3789 dmom_u_adv_u_ext[nSpace],
3790 dmom_u_adv_v_ext[nSpace],
3791 dmom_u_adv_w_ext[nSpace],
3792 mom_v_adv_ext[nSpace],
3793 dmom_v_adv_u_ext[nSpace],
3794 dmom_v_adv_v_ext[nSpace],
3795 dmom_v_adv_w_ext[nSpace],
3796 mom_w_adv_ext[nSpace],
3797 dmom_w_adv_u_ext[nSpace],
3798 dmom_w_adv_v_ext[nSpace],
3799 dmom_w_adv_w_ext[nSpace],
3800 mom_uu_diff_ten_ext[nSpace],
3801 mom_vv_diff_ten_ext[nSpace],
3802 mom_ww_diff_ten_ext[nSpace],
3803 mom_uv_diff_ten_ext[1],
3804 mom_uw_diff_ten_ext[1],
3805 mom_vu_diff_ten_ext[1],
3806 mom_vw_diff_ten_ext[1],
3807 mom_wu_diff_ten_ext[1],
3808 mom_wv_diff_ten_ext[1],
3809 mom_u_source_ext=0.0,
3810 mom_v_source_ext=0.0,
3811 mom_w_source_ext=0.0,
3813 dmom_u_ham_grad_p_ext[nSpace],
3814 dmom_u_ham_grad_u_ext[nSpace],
3816 dmom_v_ham_grad_p_ext[nSpace],
3817 dmom_v_ham_grad_v_ext[nSpace],
3819 dmom_w_ham_grad_p_ext[nSpace],
3820 dmom_w_ham_grad_w_ext[nSpace],
3821 dmom_u_adv_p_ext[nSpace],
3822 dmom_v_adv_p_ext[nSpace],
3823 dmom_w_adv_p_ext[nSpace],
3825 flux_mom_u_adv_ext=0.0,
3826 flux_mom_v_adv_ext=0.0,
3827 flux_mom_w_adv_ext=0.0,
3828 flux_mom_uu_diff_ext=0.0,
3829 flux_mom_uv_diff_ext=0.0,
3830 flux_mom_uw_diff_ext=0.0,
3831 flux_mom_vu_diff_ext=0.0,
3832 flux_mom_vv_diff_ext=0.0,
3833 flux_mom_vw_diff_ext=0.0,
3834 flux_mom_wu_diff_ext=0.0,
3835 flux_mom_wv_diff_ext=0.0,
3836 flux_mom_ww_diff_ext=0.0,
3841 bc_mom_u_acc_ext=0.0,
3842 bc_dmom_u_acc_u_ext=0.0,
3843 bc_mom_v_acc_ext=0.0,
3844 bc_dmom_v_acc_v_ext=0.0,
3845 bc_mom_w_acc_ext=0.0,
3846 bc_dmom_w_acc_w_ext=0.0,
3847 bc_mass_adv_ext[nSpace],
3848 bc_dmass_adv_u_ext[nSpace],
3849 bc_dmass_adv_v_ext[nSpace],
3850 bc_dmass_adv_w_ext[nSpace],
3851 bc_mom_u_adv_ext[nSpace],
3852 bc_dmom_u_adv_u_ext[nSpace],
3853 bc_dmom_u_adv_v_ext[nSpace],
3854 bc_dmom_u_adv_w_ext[nSpace],
3855 bc_mom_v_adv_ext[nSpace],
3856 bc_dmom_v_adv_u_ext[nSpace],
3857 bc_dmom_v_adv_v_ext[nSpace],
3858 bc_dmom_v_adv_w_ext[nSpace],
3859 bc_mom_w_adv_ext[nSpace],
3860 bc_dmom_w_adv_u_ext[nSpace],
3861 bc_dmom_w_adv_v_ext[nSpace],
3862 bc_dmom_w_adv_w_ext[nSpace],
3863 bc_mom_uu_diff_ten_ext[nSpace],
3864 bc_mom_vv_diff_ten_ext[nSpace],
3865 bc_mom_ww_diff_ten_ext[nSpace],
3866 bc_mom_uv_diff_ten_ext[1],
3867 bc_mom_uw_diff_ten_ext[1],
3868 bc_mom_vu_diff_ten_ext[1],
3869 bc_mom_vw_diff_ten_ext[1],
3870 bc_mom_wu_diff_ten_ext[1],
3871 bc_mom_wv_diff_ten_ext[1],
3872 bc_mom_u_source_ext=0.0,
3873 bc_mom_v_source_ext=0.0,
3874 bc_mom_w_source_ext=0.0,
3875 bc_mom_u_ham_ext=0.0,
3876 bc_dmom_u_ham_grad_p_ext[nSpace],
3877 bc_dmom_u_ham_grad_u_ext[nSpace],
3878 bc_mom_v_ham_ext=0.0,
3879 bc_dmom_v_ham_grad_p_ext[nSpace],
3880 bc_dmom_v_ham_grad_v_ext[nSpace],
3881 bc_mom_w_ham_ext=0.0,
3882 bc_dmom_w_ham_grad_p_ext[nSpace],
3883 bc_dmom_w_ham_grad_w_ext[nSpace],
3884 jac_ext[nSpace*nSpace],
3886 jacInv_ext[nSpace*nSpace],
3887 boundaryJac[nSpace*(nSpace-1)],
3888 metricTensor[(nSpace-1)*(nSpace-1)],
3889 metricTensorDetSqrt,
3890 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
3891 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
3892 vel_grad_test_dS[nDOF_trial_element*nSpace],
3893 normal[2],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3897 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3898 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;
3900 ck.calculateMapping_elementBoundary(eN,
3906 mesh_trial_trace_ref.data(),
3907 mesh_grad_trial_trace_ref.data(),
3908 boundaryJac_ref.data(),
3914 metricTensorDetSqrt,
3918 ck.calculateMappingVelocity_elementBoundary(eN,
3922 mesh_velocity_dof.data(),
3924 mesh_trial_trace_ref.data(),
3925 xt_ext,yt_ext,zt_ext,
3937 dS = metricTensorDetSqrt*dS_ref[kb];
3940 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3941 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
3943 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3944 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3945 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3950 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3954 p_ext = ebqe_p[ebNE_kb];
3955 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3956 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3959 for (
int I=0;I<nSpace;I++)
3960 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace + I];
3961 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3962 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3965 for (
int j=0;j<nDOF_trial_element;j++)
3968 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3969 for (
int I=0;I<nSpace;I++)
3970 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3972 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
3974 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
3975 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
3978 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
3982 double distance_to_omega_solid = 1e10;
3983 if (use_ball_as_particle == 1)
3985 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
3989 distance_to_omega_solid = ebq_global_phi_solid[ebN*nQuadraturePoints_elementBoundary+kb];
3991 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);
4000 elementDiameter[eN],
4001 smagorinskyConstant,
4002 turbulenceClosureModel,
4005 ebqe_vf_ext[ebNE_kb],
4006 ebqe_phi_ext[ebNE_kb],
4007 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
4008 distance_to_omega_solid,
4009 ebqe_kappa_phi_ext[ebNE_kb],
4021 ebqe_velocity_star[ebNE_kb_nSpace+0],
4022 ebqe_velocity_star[ebNE_kb_nSpace+1],
4023 ebqe_velocity_star[ebNE_kb_nSpace+1],
4047 mom_uu_diff_ten_ext,
4048 mom_vv_diff_ten_ext,
4049 mom_ww_diff_ten_ext,
4050 mom_uv_diff_ten_ext,
4051 mom_uw_diff_ten_ext,
4052 mom_vu_diff_ten_ext,
4053 mom_vw_diff_ten_ext,
4054 mom_wu_diff_ten_ext,
4055 mom_wv_diff_ten_ext,
4060 dmom_u_ham_grad_p_ext,
4061 dmom_u_ham_grad_u_ext,
4063 dmom_v_ham_grad_p_ext,
4064 dmom_v_ham_grad_v_ext,
4066 dmom_w_ham_grad_p_ext,
4067 dmom_w_ham_grad_w_ext,
4075 MATERIAL_PARAMETERS_AS_FUNCTION,
4076 ebqe_density_as_function[ebNE_kb],
4077 ebqe_dynamic_viscosity_as_function[ebNE_kb],
4080 use_ball_as_particle,
4083 ball_velocity.data(),
4084 ball_angular_velocity.data(),
4085 INT_BY_PARTS_PRESSURE);
4094 elementDiameter[eN],
4095 smagorinskyConstant,
4096 turbulenceClosureModel,
4099 bc_ebqe_vf_ext[ebNE_kb],
4100 bc_ebqe_phi_ext[ebNE_kb],
4101 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
4102 distance_to_omega_solid,
4103 ebqe_kappa_phi_ext[ebNE_kb],
4115 ebqe_velocity_star[ebNE_kb_nSpace+0],
4116 ebqe_velocity_star[ebNE_kb_nSpace+1],
4117 ebqe_velocity_star[ebNE_kb_nSpace+1],
4118 bc_eddy_viscosity_ext,
4120 bc_dmom_u_acc_u_ext,
4122 bc_dmom_v_acc_v_ext,
4124 bc_dmom_w_acc_w_ext,
4130 bc_dmom_u_adv_u_ext,
4131 bc_dmom_u_adv_v_ext,
4132 bc_dmom_u_adv_w_ext,
4134 bc_dmom_v_adv_u_ext,
4135 bc_dmom_v_adv_v_ext,
4136 bc_dmom_v_adv_w_ext,
4138 bc_dmom_w_adv_u_ext,
4139 bc_dmom_w_adv_v_ext,
4140 bc_dmom_w_adv_w_ext,
4141 bc_mom_uu_diff_ten_ext,
4142 bc_mom_vv_diff_ten_ext,
4143 bc_mom_ww_diff_ten_ext,
4144 bc_mom_uv_diff_ten_ext,
4145 bc_mom_uw_diff_ten_ext,
4146 bc_mom_vu_diff_ten_ext,
4147 bc_mom_vw_diff_ten_ext,
4148 bc_mom_wu_diff_ten_ext,
4149 bc_mom_wv_diff_ten_ext,
4150 bc_mom_u_source_ext,
4151 bc_mom_v_source_ext,
4152 bc_mom_w_source_ext,
4154 bc_dmom_u_ham_grad_p_ext,
4155 bc_dmom_u_ham_grad_u_ext,
4157 bc_dmom_v_ham_grad_p_ext,
4158 bc_dmom_v_ham_grad_v_ext,
4160 bc_dmom_w_ham_grad_p_ext,
4161 bc_dmom_w_ham_grad_w_ext,
4169 MATERIAL_PARAMETERS_AS_FUNCTION,
4170 ebqe_density_as_function[ebNE_kb],
4171 ebqe_dynamic_viscosity_as_function[ebNE_kb],
4174 use_ball_as_particle,
4177 ball_velocity.data(),
4178 ball_angular_velocity.data(),
4179 INT_BY_PARTS_PRESSURE);
4182 if (turbulenceClosureModel >= 3)
4184 const double turb_var_grad_0_dummy[2] = {0.,0.};
4185 const double c_mu = 0.09;
4194 ebqe_vf_ext[ebNE_kb],
4195 ebqe_phi_ext[ebNE_kb],
4198 ebqe_turb_var_0[ebNE_kb],
4199 ebqe_turb_var_1[ebNE_kb],
4200 turb_var_grad_0_dummy,
4202 mom_uu_diff_ten_ext,
4203 mom_vv_diff_ten_ext,
4204 mom_ww_diff_ten_ext,
4205 mom_uv_diff_ten_ext,
4206 mom_uw_diff_ten_ext,
4207 mom_vu_diff_ten_ext,
4208 mom_vw_diff_ten_ext,
4209 mom_wu_diff_ten_ext,
4210 mom_wv_diff_ten_ext,
4223 bc_ebqe_vf_ext[ebNE_kb],
4224 bc_ebqe_phi_ext[ebNE_kb],
4227 ebqe_turb_var_0[ebNE_kb],
4228 ebqe_turb_var_1[ebNE_kb],
4229 turb_var_grad_0_dummy,
4230 bc_eddy_viscosity_ext,
4231 bc_mom_uu_diff_ten_ext,
4232 bc_mom_vv_diff_ten_ext,
4233 bc_mom_ww_diff_ten_ext,
4234 bc_mom_uv_diff_ten_ext,
4235 bc_mom_uw_diff_ten_ext,
4236 bc_mom_vu_diff_ten_ext,
4237 bc_mom_vw_diff_ten_ext,
4238 bc_mom_wu_diff_ten_ext,
4239 bc_mom_wv_diff_ten_ext,
4240 bc_mom_u_source_ext,
4241 bc_mom_v_source_ext,
4242 bc_mom_w_source_ext);
4249 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
4250 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
4252 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
4253 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
4256 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
4257 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
4259 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
4260 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
4272 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
4273 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
4276 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
4277 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
4286 ck.calculateGScale(G,normal,h_penalty);
4287 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
4289 isDOFBoundary_u[ebNE_kb],
4290 isDOFBoundary_v[ebNE_kb],
4291 isDOFBoundary_w[ebNE_kb],
4292 isAdvectiveFluxBoundary_p[ebNE_kb],
4293 isAdvectiveFluxBoundary_u[ebNE_kb],
4294 isAdvectiveFluxBoundary_v[ebNE_kb],
4295 isAdvectiveFluxBoundary_w[ebNE_kb],
4296 dmom_u_ham_grad_p_ext[0],
4297 bc_dmom_u_ham_grad_p_ext[0],
4308 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),
4309 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
4310 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
4311 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
4339 &ebqe_velocity_star[ebNE_kb_nSpace],
4340 &ebqe_velocity[ebNE_kb_nSpace]);
4342 for (
int I=0;I<nSpace;I++)
4344 ebqe_grad_u[ebNE_kb_nSpace+I] = grad_u_ext[I];
4345 ebqe_grad_v[ebNE_kb_nSpace+I] = grad_v_ext[I];
4349 ebqe_phi_ext[ebNE_kb],
4350 sdInfo_u_u_rowptr.data(),
4351 sdInfo_u_u_colind.data(),
4352 isDOFBoundary_u[ebNE_kb],
4353 isDiffusiveFluxBoundary_u[ebNE_kb],
4355 bc_mom_uu_diff_ten_ext,
4357 ebqe_bc_flux_u_diff_ext[ebNE_kb],
4358 mom_uu_diff_ten_ext,
4362 flux_mom_uu_diff_ext);
4364 ebqe_phi_ext[ebNE_kb],
4365 sdInfo_u_v_rowptr.data(),
4366 sdInfo_u_v_colind.data(),
4367 isDOFBoundary_v[ebNE_kb],
4368 isDiffusiveFluxBoundary_v[ebNE_kb],
4370 bc_mom_uv_diff_ten_ext,
4373 mom_uv_diff_ten_ext,
4377 flux_mom_uv_diff_ext);
4394 ebqe_phi_ext[ebNE_kb],
4395 sdInfo_v_u_rowptr.data(),
4396 sdInfo_v_u_colind.data(),
4397 isDOFBoundary_u[ebNE_kb],
4398 isDiffusiveFluxBoundary_u[ebNE_kb],
4400 bc_mom_vu_diff_ten_ext,
4403 mom_vu_diff_ten_ext,
4407 flux_mom_vu_diff_ext);
4409 ebqe_phi_ext[ebNE_kb],
4410 sdInfo_v_v_rowptr.data(),
4411 sdInfo_v_v_colind.data(),
4412 isDOFBoundary_v[ebNE_kb],
4413 isDiffusiveFluxBoundary_v[ebNE_kb],
4415 bc_mom_vv_diff_ten_ext,
4417 ebqe_bc_flux_v_diff_ext[ebNE_kb],
4418 mom_vv_diff_ten_ext,
4422 flux_mom_vv_diff_ext);
4483 flux[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
4491 if (ebN < nElementBoundaries_owned)
4493 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];
4494 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];
4497 force_p_x = p_ext*normal[0];
4498 force_p_y = p_ext*normal[1];
4501 force_x = force_p_x + force_v_x;
4502 force_y = force_p_y + force_v_y;
4505 r_x = x_ext - barycenters[3*boundaryFlags[ebN]+0];
4506 r_y = y_ext - barycenters[3*boundaryFlags[ebN]+1];
4509 wettedAreas[boundaryFlags[ebN]] += dS*(1.0-ebqe_vf_ext[ebNE_kb]);
4511 netForces_p[3*boundaryFlags[ebN]+0] += force_p_x*dS;
4512 netForces_p[3*boundaryFlags[ebN]+1] += force_p_y*dS;
4515 netForces_v[3*boundaryFlags[ebN]+0] += force_v_x*dS;
4516 netForces_v[3*boundaryFlags[ebN]+1] += force_v_y*dS;
4521 netMoments[3*boundaryFlags[ebN]+2] += (r_x*force_y - r_y*force_x)*dS;
4526 for (
int i=0;i<nDOF_test_element;i++)
4532 elementResidual_u[i] +=
4533 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[0] : 0.) +
4534 ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
4535 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
4536 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
4537 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4538 isDiffusiveFluxBoundary_u[ebNE_kb],
4543 sdInfo_u_u_rowptr.data(),
4544 sdInfo_u_u_colind.data(),
4545 mom_uu_diff_ten_ext,
4546 &vel_grad_test_dS[i*nSpace])+
4547 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4548 isDiffusiveFluxBoundary_u[ebNE_kb],
4553 sdInfo_u_v_rowptr.data(),
4554 sdInfo_u_v_colind.data(),
4555 mom_uv_diff_ten_ext,
4556 &vel_grad_test_dS[i*nSpace]);
4567 elementResidual_v[i] +=
4568 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[1] : 0.) +
4569 ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
4570 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
4571 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
4572 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4573 isDiffusiveFluxBoundary_v[ebNE_kb],
4578 sdInfo_v_u_rowptr.data(),
4579 sdInfo_v_u_colind.data(),
4580 mom_vu_diff_ten_ext,
4581 &vel_grad_test_dS[i*nSpace])+
4582 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4583 isDiffusiveFluxBoundary_v[ebNE_kb],
4588 sdInfo_v_v_rowptr.data(),
4589 sdInfo_v_v_colind.data(),
4590 mom_vv_diff_ten_ext,
4591 &vel_grad_test_dS[i*nSpace]);
4644 for (
int i=0;i<nDOF_test_element;i++)
4646 int eN_i = eN*nDOF_test_element+i;
4651 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=elementResidual_u[i];
4652 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=elementResidual_v[i];
4664 particle_surfaceArea[0] = cut_cell_boundary_length;
4665 particle_netForces[(0+nParticles)*3 +0] = p_force_x;
4666 particle_netForces[(0+nParticles)*3 +1] = p_force_y;
4667 std::cout<<
"===end mesh==="<<std::endl<<std::flush;
4674 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4675 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4676 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4677 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
4678 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
4679 double PSTAB = args.
scalar<
double>(
"PSTAB");
4680 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4681 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
4682 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4683 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4684 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4685 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4686 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4687 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
4688 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
4689 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
4690 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
4691 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4692 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4693 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
4694 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4695 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4696 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4697 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4698 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4699 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4700 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4701 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4702 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4703 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4704 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4705 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4706 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4707 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4708 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4709 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
4710 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4711 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4712 double hFactor = args.
scalar<
double>(
"hFactor");
4713 int nElements_global = args.
scalar<
int>(
"nElements_global");
4714 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
4715 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
4716 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
4717 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
4718 double useRBLES = args.
scalar<
double>(
"useRBLES");
4719 double useMetrics = args.
scalar<
double>(
"useMetrics");
4720 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4721 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4722 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4723 double sigma = args.
scalar<
double>(
"sigma");
4728 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4729 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4730 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4731 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4732 double C_dg = args.
scalar<
double>(
"C_dg");
4733 double C_b = args.
scalar<
double>(
"C_b");
4734 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
4735 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
4736 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
4737 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
4738 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
4739 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
4740 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
4741 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
4742 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
4743 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
4744 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
4745 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
4746 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
4747 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
4748 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
4749 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
4750 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
4751 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4752 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4753 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4754 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4755 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4756 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4757 xt::pyarray<double>& g = args.
array<
double>(
"g");
4758 const double useVF = args.
scalar<
double>(
"useVF");
4759 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4760 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4761 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
4762 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4763 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4764 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4765 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4766 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4767 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
4768 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
4769 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4770 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
4771 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4772 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4773 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4774 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4775 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4776 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4777 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4778 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4779 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4780 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4781 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4782 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4783 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4784 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4785 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4786 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4787 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4788 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4789 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4790 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4791 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4792 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4793 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
4794 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
4795 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
4796 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
4797 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
4798 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
4799 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
4800 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
4801 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
4802 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
4803 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
4804 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
4805 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
4806 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
4807 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
4808 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
4809 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
4810 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
4811 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
4812 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
4813 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
4814 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
4815 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
4816 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
4817 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
4818 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
4819 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
4820 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
4821 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
4822 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
4823 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
4824 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
4825 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
4826 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4827 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4828 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
4829 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4830 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4831 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4832 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
4833 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4834 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
4835 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4836 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4837 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
4838 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
4839 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
4840 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4841 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4842 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4843 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4844 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4845 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4846 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4847 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4848 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4849 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4850 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4851 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4852 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4853 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4854 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4855 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4856 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4857 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4858 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4859 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4860 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4861 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4862 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4863 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
4864 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
4865 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
4866 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
4867 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
4868 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
4869 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
4870 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
4871 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
4872 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
4873 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
4874 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
4875 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
4876 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
4877 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
4878 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
4879 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
4880 int nParticles = args.
scalar<
int>(
"nParticles");
4881 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
4882 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
4883 double particle_beta = args.
scalar<
double>(
"particle_beta");
4884 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
4885 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
4886 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
4887 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
4888 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
4889 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
4890 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
4891 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
4892 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
4893 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
4894 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
4895 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
4896 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
4897 double dt = args.
scalar<
double>(
"dt");
4898 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
4899 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
4900 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
4901 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
4902 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
4903 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
4904 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
4905 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
4906 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
4907 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
4908 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
4909 int offset_u = args.
scalar<
int>(
"offset_u");
4910 int offset_v = args.
scalar<
int>(
"offset_v");
4911 int offset_w = args.
scalar<
int>(
"offset_w");
4912 int stride_u = args.
scalar<
int>(
"stride_u");
4913 int stride_v = args.
scalar<
int>(
"stride_v");
4914 int stride_w = args.
scalar<
int>(
"stride_w");
4915 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
4916 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
4917 xt::pyarray<int>& rowptr = args.
array<
int>(
"rowptr");
4918 xt::pyarray<int>& colind = args.
array<
int>(
"colind");
4919 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
4923 std::valarray<double> particle_surfaceArea(nParticles), particle_netForces(nParticles*3*3), particle_netMoments(nParticles*3);
4924 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
4925 for(
int eN=0;eN<nElements_global;eN++)
4927 double eps_rho,eps_mu;
4928 double element_active=1.0;
4930 double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
4931 elementJacobian_p_u[nDOF_test_element][nDOF_trial_element],
4932 elementJacobian_p_v[nDOF_test_element][nDOF_trial_element],
4933 elementJacobian_p_w[nDOF_test_element][nDOF_trial_element],
4934 elementJacobian_u_p[nDOF_test_element][nDOF_trial_element],
4935 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
4936 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
4937 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
4938 elementJacobian_v_p[nDOF_test_element][nDOF_trial_element],
4939 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
4940 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
4941 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
4942 elementJacobian_w_p[nDOF_test_element][nDOF_trial_element],
4943 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
4944 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
4945 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
4946 for (
int i=0;i<nDOF_test_element;i++)
4947 for (
int j=0;j<nDOF_trial_element;j++)
4949 elementJacobian_p_p[i][j]=0.0;
4950 elementJacobian_p_u[i][j]=0.0;
4951 elementJacobian_p_v[i][j]=0.0;
4952 elementJacobian_p_w[i][j]=0.0;
4953 elementJacobian_u_p[i][j]=0.0;
4954 elementJacobian_u_u[i][j]=0.0;
4955 elementJacobian_u_v[i][j]=0.0;
4956 elementJacobian_u_w[i][j]=0.0;
4957 elementJacobian_v_p[i][j]=0.0;
4958 elementJacobian_v_u[i][j]=0.0;
4959 elementJacobian_v_v[i][j]=0.0;
4960 elementJacobian_v_w[i][j]=0.0;
4961 elementJacobian_w_p[i][j]=0.0;
4962 elementJacobian_w_u[i][j]=0.0;
4963 elementJacobian_w_v[i][j]=0.0;
4964 elementJacobian_w_w[i][j]=0.0;
4975 double _distance[nDOF_mesh_trial_element]={0.0};
4977 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4979 if(use_ball_as_particle==1)
4982 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
4983 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
4984 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
4989 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
4991 if ( _distance[I] >= 0)
4994 if (pos_counter == 2)
4999 for (
int I=0;I<nDOF_mesh_trial_element;I++)
5001 if (_distance[I] < 0)
5004 assert(opp_node >=0);
5005 assert(opp_node <nDOF_mesh_trial_element);
5007 else if (pos_counter == 3)
5016 if(use_ball_as_particle==1)
5018 for (
int I=0;I<nDOF_mesh_trial_element;I++)
5020 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
5021 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
5022 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
5023 phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]]);
5025 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
5026 for (
int j=0;j<nDOF_mesh_trial_element;j++)
5028 int eN_j = eN*nDOF_mesh_trial_element+j;
5029 element_phi[j] = phi_dof[p_l2g[eN_j]];
5030 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
5032 double element_nodes[nDOF_mesh_trial_element*3];
5033 for (
int i=0;i<nDOF_mesh_trial_element;i++)
5035 int eN_i=eN*nDOF_mesh_trial_element+i;
5036 for(
int I=0;I<3;I++)
5037 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
5039 gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
5040 gf.
calculate(element_phi, element_nodes, x_ref.data(),
false);
5041 for (
int k=0;k<nQuadraturePoints_element;k++)
5045 int eN_k = eN*nQuadraturePoints_element+k,
5046 eN_k_nSpace = eN_k*nSpace,
5048 eN_nDOF_trial_element = eN*nDOF_trial_element;
5051 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
5052 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
5061 dmass_adv_u[nSpace],
5062 dmass_adv_v[nSpace],
5063 dmass_adv_w[nSpace],
5065 dmom_u_adv_u[nSpace],
5066 dmom_u_adv_v[nSpace],
5067 dmom_u_adv_w[nSpace],
5069 dmom_v_adv_u[nSpace],
5070 dmom_v_adv_v[nSpace],
5071 dmom_v_adv_w[nSpace],
5073 dmom_w_adv_u[nSpace],
5074 dmom_w_adv_v[nSpace],
5075 dmom_w_adv_w[nSpace],
5076 mom_uu_diff_ten[nSpace],
5077 mom_vv_diff_ten[nSpace],
5078 mom_ww_diff_ten[nSpace],
5089 dmom_u_ham_grad_p[nSpace],
5090 dmom_u_ham_grad_u[nSpace],
5092 dmom_v_ham_grad_p[nSpace],
5093 dmom_v_ham_grad_v[nSpace],
5095 dmom_w_ham_grad_p[nSpace],
5096 dmom_w_ham_grad_w[nSpace],
5107 dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
5108 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
5109 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
5110 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
5111 Lstar_u_p[nDOF_test_element],
5112 Lstar_v_p[nDOF_test_element],
5113 Lstar_w_p[nDOF_test_element],
5114 Lstar_u_u[nDOF_test_element],
5115 Lstar_v_v[nDOF_test_element],
5116 Lstar_w_w[nDOF_test_element],
5117 Lstar_p_u[nDOF_test_element],
5118 Lstar_p_v[nDOF_test_element],
5119 Lstar_p_w[nDOF_test_element],
5124 dsubgridError_p_u[nDOF_trial_element],
5125 dsubgridError_p_v[nDOF_trial_element],
5126 dsubgridError_p_w[nDOF_trial_element],
5127 dsubgridError_u_p[nDOF_trial_element],
5128 dsubgridError_u_u[nDOF_trial_element],
5129 dsubgridError_v_p[nDOF_trial_element],
5130 dsubgridError_v_v[nDOF_trial_element],
5131 dsubgridError_w_p[nDOF_trial_element],
5132 dsubgridError_w_w[nDOF_trial_element],
5133 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
5134 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
5137 jacInv[nSpace*nSpace],
5138 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
5139 vel_hess_trial[nDOF_trial_element*
nSpace2],
5141 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_test_element],
5142 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
5147 dmom_u_source[nSpace],
5148 dmom_v_source[nSpace],
5149 dmom_w_source[nSpace],
5152 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
5154 ck.calculateMapping_element(eN,
5158 mesh_trial_ref.data(),
5159 mesh_grad_trial_ref.data(),
5164 ck.calculateH_element(eN,
5166 nodeDiametersArray.data(),
5168 mesh_trial_ref.data(),
5170 ck.calculateMappingVelocity_element(eN,
5172 mesh_velocity_dof.data(),
5174 mesh_trial_ref.data(),
5179 dV = fabs(jacDet)*dV_ref[k];
5180 ck.calculateG(jacInv,G,G_dd_G,tr_G);
5183 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
5184 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
5185 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
5189 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
5190 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
5194 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
5195 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
5199 for (
int I=0;I<nSpace;I++)
5200 grad_p[I] = q_grad_p[eN_k_nSpace+I];
5201 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
5202 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
5203 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
5204 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
5207 for (
int j=0;j<nDOF_trial_element;j++)
5210 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
5211 for (
int I=0;I<nSpace;I++)
5214 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
5218 double div_mesh_velocity=0.0;
5219 int NDOF_MESH_TRIAL_ELEMENT=3;
5220 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
5222 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
5223 div_mesh_velocity +=
5224 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*2+0] +
5225 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*2+1];
5227 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
5230 porosity = 1.0 - q_vos[eN_k];
5235 double distance_to_omega_solid = phi_solid[eN_k];
5236 double eddy_viscosity(0.),rhoSave,nuSave;
5245 elementDiameter[eN],
5246 smagorinskyConstant,
5247 turbulenceClosureModel,
5252 &normal_phi[eN_k_nSpace],
5253 distance_to_omega_solid,
5266 q_velocity_sge[eN_k_nSpace+0],
5267 q_velocity_sge[eN_k_nSpace+1],
5268 q_velocity_sge[eN_k_nSpace+1],
5320 MATERIAL_PARAMETERS_AS_FUNCTION,
5321 density_as_function[eN_k],
5322 dynamic_viscosity_as_function[eN_k],
5325 use_ball_as_particle,
5328 ball_velocity.data(),
5329 ball_angular_velocity.data(),
5330 INT_BY_PARTS_PRESSURE);
5332 mass_source = q_mass_source[eN_k];
5333 for (
int I=0;I<nSpace;I++)
5335 dmom_u_source[I] = 0.0;
5336 dmom_v_source[I] = 0.0;
5337 dmom_w_source[I] = 0.0;
5358 q_velocity_sge[eN_k_nSpace+0],
5359 q_velocity_sge[eN_k_nSpace+1],
5360 q_velocity_sge[eN_k_nSpace+1],
5361 eps_solid[elementFlags[eN]],
5363 q_velocity_solid[eN_k_nSpace+0],
5364 q_velocity_solid[eN_k_nSpace+1],
5365 q_velocity_solid[eN_k_nSpace+1],
5366 q_velocityStar_solid[eN_k_nSpace+0],
5367 q_velocityStar_solid[eN_k_nSpace+1],
5368 q_velocityStar_solid[eN_k_nSpace+1],
5375 q_grad_vos[eN_k_nSpace+0],
5376 q_grad_vos[eN_k_nSpace+1],
5377 q_grad_vos[eN_k_nSpace+1]);
5379 double C_particles=0.0;
5380 if(nParticles > 0 && USE_SBM==0)
5385 nQuadraturePoints_global,
5386 &particle_signed_distances[eN_k],
5387 &particle_signed_distance_normals[eN_k_3d],
5388 &particle_velocities[eN_k_3d],
5389 particle_centroids.data(),
5390 use_ball_as_particle,
5393 ball_velocity.data(),
5394 ball_angular_velocity.data(),
5396 particle_penalty_constant/h_phi,
5397 particle_alpha/h_phi,
5398 particle_beta/h_phi,
5415 q_velocity_sge[eN_k_nSpace+0],
5416 q_velocity_sge[eN_k_nSpace+1],
5417 q_velocity_sge[eN_k_nSpace+1],
5440 &particle_netForces[0],
5441 &particle_netMoments[0],
5442 &particle_surfaceArea[0]);
5444 if (turbulenceClosureModel >= 3)
5446 const double c_mu = 0.09;
5461 &q_turb_var_grad_0[eN_k_nSpace],
5481 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
5482 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
5484 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
5485 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
5488 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
5489 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
5491 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
5492 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
5505 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5511 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5525 mom_u_acc_t *= dmom_u_acc_u;
5526 mom_v_acc_t *= dmom_v_acc_v;
5529 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
5530 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
5536 ck.Mass_strong(-q_dvos_dt[eN_k]) +
5537 ck.Advection_strong(dmass_adv_u,grad_u) +
5538 ck.Advection_strong(dmass_adv_v,grad_v) +
5540 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
5542 ck.Reaction_strong(mass_source);
5546 ck.Mass_strong(mom_u_acc_t) +
5547 ck.Advection_strong(dmom_adv_sge,grad_u) +
5548 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
5549 ck.Reaction_strong(mom_u_source) -
5550 ck.Reaction_strong(
u*div_mesh_velocity);
5553 ck.Mass_strong(mom_v_acc_t) +
5554 ck.Advection_strong(dmom_adv_sge,grad_v) +
5555 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
5556 ck.Reaction_strong(mom_v_source) -
5557 ck.Reaction_strong(
v*div_mesh_velocity);
5567 for (
int j=0;j<nDOF_trial_element;j++)
5569 int j_nSpace = j*nSpace;
5570 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
5571 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
5574 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
5575 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5576 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5577 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5579 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
5580 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5581 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5582 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5590 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j]);
5591 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j]);
5597 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
5599 elementDiameter[eN],
5604 dmom_u_ham_grad_p[0],
5614 dmom_u_ham_grad_p[0],
5620 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
5621 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
5654 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
5655 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
5659 for (
int i=0;i<nDOF_test_element;i++)
5661 int i_nSpace = i*nSpace;
5662 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
5663 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
5665 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5666 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5668 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5669 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5672 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
5673 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
5678 dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5679 dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5682 dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5683 dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5691 double unit_normal[nSpace];
5692 double norm_grad_phi = 0.;
5693 for (
int I=0;I<nSpace;I++)
5694 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
5695 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
5696 for (
int I=0;I<nSpace;I++)
5697 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
5698 double delta =
gf.
D(eps_mu,
phi[eN_k]);
5699 double vel_tgrad_test_i[nSpace], vel_tgrad_test_j[nSpace];
5703 for(
int i=0;i<nDOF_test_element;i++)
5705 int i_nSpace = i*nSpace;
5707 &vel_grad_trial[i_nSpace],
5709 for(
int j=0;j<nDOF_trial_element;j++)
5711 int j_nSpace = j*nSpace;
5713 &vel_grad_trial[j_nSpace],
5729 elementJacobian_u_u[i][j] +=
5730 ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5731 ck.HamiltonianJacobian_weak(dmom_u_ham_grad_u,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5732 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5733 ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5735 ck.ReactionJacobian_weak(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5738 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
5739 ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5741 ck.NumericalDiffusion(dt*delta*sigma*dV,
5745 elementJacobian_u_v[i][j] +=
5746 ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5747 ck.SimpleDiffusionJacobian_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5749 ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5761 elementJacobian_v_u[i][j] +=
5762 ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5763 ck.SimpleDiffusionJacobian_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5765 ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5768 elementJacobian_v_v[i][j] +=
5769 ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5770 ck.HamiltonianJacobian_weak(dmom_v_ham_grad_v,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5771 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5772 ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5774 ck.ReactionJacobian_weak(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5777 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
5778 ck.NumericalDiffusionJacobian(q_numDiff_v_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5780 ck.NumericalDiffusion(dt*delta*sigma*dV,
5821 for (
int i=0;i<nDOF_test_element;i++)
5823 int eN_i = eN*nDOF_test_element+i;
5824 for (
int j=0;j<nDOF_trial_element;j++)
5826 int eN_i_j = eN_i*nDOF_trial_element+j;
5833 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += element_active*elementJacobian_u_u[i][j];
5834 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_u_v[eN_i_j]] += element_active*elementJacobian_u_v[i][j];
5838 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_v_u[eN_i_j]] += element_active*elementJacobian_v_u[i][j];
5839 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_v_v[eN_i_j]] += element_active*elementJacobian_v_v[i][j];
5851 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
5854 for (
int i=0; i<numDOFs_1D; i++)
5857 int u_gi = offset_u+stride_u*i;
5858 int v_gi = offset_v+stride_v*i;
5861 int u_ith_row_ptr = rowptr[u_gi];
5862 int v_ith_row_ptr = rowptr[v_gi];
5865 int numDOFs_ith_row = rowptr_1D[i+1]-rowptr_1D[i];
5866 for (
int counter = 0; counter < numDOFs_ith_row; counter++)
5869 int uu_ij = u_ith_row_ptr + (offset_u + counter*stride_u);
5870 int vv_ij = v_ith_row_ptr + (offset_v + counter*stride_v);
5873 double uStar_dij = uStar_dMatrix[ij];
5874 double vStar_dij = vStar_dMatrix[ij];
5877 globalJacobian[uu_ij] -= uStar_dij;
5878 globalJacobian[vv_ij] -= vStar_dij;
5895 eN_nDOF_trial_element = eN*nDOF_trial_element;
5896 double eps_rho,eps_mu;
5900 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5902 int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
5903 ebN_kb_nSpace = ebN_kb*nSpace,
5904 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5905 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5913 jac_ext[nSpace*nSpace],
5915 jacInv_ext[nSpace*nSpace],
5916 boundaryJac[nSpace*(nSpace-1)],
5917 metricTensor[(nSpace-1)*(nSpace-1)],
5918 metricTensorDetSqrt,
5919 vel_grad_trial_trace[nDOF_trial_element*nSpace],
5921 vel_test_dS[nDOF_test_element],
5923 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
5924 vel_grad_test_dS[nDOF_trial_element*nSpace],
5925 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
5926 ck.calculateMapping_elementBoundary(eN,
5932 mesh_trial_trace_ref.data(),
5933 mesh_grad_trial_trace_ref.data(),
5934 boundaryJac_ref.data(),
5940 metricTensorDetSqrt,
5944 ck.calculateMappingVelocity_elementBoundary(eN,
5948 mesh_velocity_dof.data(),
5950 mesh_trial_trace_ref.data(),
5951 xt_ext,yt_ext,zt_ext,
5956 dS = metricTensorDetSqrt*dS_ref[kb];
5957 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
5960 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
5962 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
5963 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
5965 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
5966 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
5968 for (
int j=0;j<nDOF_trial_element;j++)
5970 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
5971 for (
int I=0;I<nSpace;I++)
5972 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
5979 ck.calculateGScale(G,normal,h_penalty);
5985 double distance[2], P_normal[2], P_tangent[2];
5987 if(use_ball_as_particle==1)
5995 P_normal[0],P_normal[1]);
5997 ball_velocity.data(), ball_angular_velocity.data(),
5999 x_ext-dist*P_normal[0],
6000 y_ext-dist*P_normal[1],
6006 dist = ebq_global_phi_solid[ebN_kb];
6007 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*3+0];
6008 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*3+1];
6009 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*3+0];
6010 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*3+1];
6012 distance[0] = -P_normal[0]*dist;
6013 distance[1] = -P_normal[1]*dist;
6014 P_tangent[0]= -P_normal[1];
6015 P_tangent[1]= P_normal[0];
6016 assert(h_penalty>0.0);
6017 if (h_penalty < std::abs(dist))
6018 h_penalty = std::abs(dist);
6021 double C_adim =
C_sbm*visco/h_penalty;
6022 double beta_adim =
beta_sbm*visco/h_penalty;
6024 for (
int i=0;i<nDOF_test_element;i++)
6026 int eN_i = eN*nDOF_test_element+i;
6027 double phi_i = vel_test_dS[i];
6028 double* grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
6030 const double grad_phi_i_dot_t =
get_dot_product(P_tangent,grad_phi_i);
6033 const double zero_vec[2]={0.,0.};
6034 for (
int j=0;j<nDOF_trial_element;j++)
6039 + i*nDOF_trial_element
6042 double phi_j = vel_test_dS[j]/dS;
6043 const double grad_phi_j[2]={vel_grad_test_dS[j*nSpace+0]/dS,
6044 vel_grad_test_dS[j*nSpace+1]/dS};
6045 const double grad_phi_j_dot_d =
get_dot_product(distance, grad_phi_j);
6046 const double grad_phi_j_dot_t =
get_dot_product(P_tangent,grad_phi_j);
6050 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6052 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6057 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
6058 visco * phi_i * res[0];
6059 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
6060 visco * phi_i * res[1];
6063 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
6064 visco * phi_i * res[0];
6065 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
6066 visco * phi_i * res[1];
6070 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
6071 visco * phi_j * res[0];
6072 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
6073 visco * phi_j * res[1];
6075 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
6076 visco * phi_j * res[0];
6077 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
6078 visco * phi_j * res[1];
6081 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6082 C_adim*grad_phi_i_dot_d*phi_j;
6083 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6084 C_adim*grad_phi_i_dot_d*phi_j;
6087 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6088 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
6089 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6090 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
6093 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6094 C_adim*grad_phi_j_dot_d*phi_i;
6095 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6096 C_adim*grad_phi_j_dot_d*phi_i;
6100 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
6101 visco * grad_phi_j_dot_d * res[0];
6102 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
6103 visco * grad_phi_j_dot_d * res[1];
6106 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
6107 visco * grad_phi_j_dot_d * res[0] ;
6108 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
6109 visco * grad_phi_j_dot_d * res[1];
6114 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
6115 beta_adim*grad_phi_j_dot_t*grad_phi_i_dot_t;
6116 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
6117 beta_adim*grad_phi_j_dot_t*grad_phi_i_dot_t;
6130 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
6132 int ebN = exteriorElementBoundariesArray[ebNE],
6133 eN = elementBoundaryElementsArray[ebN*2+0],
6134 eN_nDOF_trial_element = eN*nDOF_trial_element,
6135 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0];
6136 double eps_rho,eps_mu;
6137 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
6139 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
6140 ebNE_kb_nSpace = ebNE_kb*nSpace,
6141 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
6142 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
6153 dmom_u_acc_u_ext=0.0,
6155 dmom_v_acc_v_ext=0.0,
6157 dmom_w_acc_w_ext=0.0,
6158 mass_adv_ext[nSpace],
6159 dmass_adv_u_ext[nSpace],
6160 dmass_adv_v_ext[nSpace],
6161 dmass_adv_w_ext[nSpace],
6162 mom_u_adv_ext[nSpace],
6163 dmom_u_adv_u_ext[nSpace],
6164 dmom_u_adv_v_ext[nSpace],
6165 dmom_u_adv_w_ext[nSpace],
6166 mom_v_adv_ext[nSpace],
6167 dmom_v_adv_u_ext[nSpace],
6168 dmom_v_adv_v_ext[nSpace],
6169 dmom_v_adv_w_ext[nSpace],
6170 mom_w_adv_ext[nSpace],
6171 dmom_w_adv_u_ext[nSpace],
6172 dmom_w_adv_v_ext[nSpace],
6173 dmom_w_adv_w_ext[nSpace],
6174 mom_uu_diff_ten_ext[nSpace],
6175 mom_vv_diff_ten_ext[nSpace],
6176 mom_ww_diff_ten_ext[nSpace],
6177 mom_uv_diff_ten_ext[1],
6178 mom_uw_diff_ten_ext[1],
6179 mom_vu_diff_ten_ext[1],
6180 mom_vw_diff_ten_ext[1],
6181 mom_wu_diff_ten_ext[1],
6182 mom_wv_diff_ten_ext[1],
6183 mom_u_source_ext=0.0,
6184 mom_v_source_ext=0.0,
6185 mom_w_source_ext=0.0,
6187 dmom_u_ham_grad_p_ext[nSpace],
6188 dmom_u_ham_grad_u_ext[nSpace],
6190 dmom_v_ham_grad_p_ext[nSpace],
6191 dmom_v_ham_grad_v_ext[nSpace],
6193 dmom_w_ham_grad_p_ext[nSpace],
6194 dmom_w_ham_grad_w_ext[nSpace],
6195 dmom_u_adv_p_ext[nSpace],
6196 dmom_v_adv_p_ext[nSpace],
6197 dmom_w_adv_p_ext[nSpace],
6198 dflux_mass_u_ext=0.0,
6199 dflux_mass_v_ext=0.0,
6200 dflux_mass_w_ext=0.0,
6201 dflux_mom_u_adv_p_ext=0.0,
6202 dflux_mom_u_adv_u_ext=0.0,
6203 dflux_mom_u_adv_v_ext=0.0,
6204 dflux_mom_u_adv_w_ext=0.0,
6205 dflux_mom_v_adv_p_ext=0.0,
6206 dflux_mom_v_adv_u_ext=0.0,
6207 dflux_mom_v_adv_v_ext=0.0,
6208 dflux_mom_v_adv_w_ext=0.0,
6209 dflux_mom_w_adv_p_ext=0.0,
6210 dflux_mom_w_adv_u_ext=0.0,
6211 dflux_mom_w_adv_v_ext=0.0,
6212 dflux_mom_w_adv_w_ext=0.0,
6217 bc_mom_u_acc_ext=0.0,
6218 bc_dmom_u_acc_u_ext=0.0,
6219 bc_mom_v_acc_ext=0.0,
6220 bc_dmom_v_acc_v_ext=0.0,
6221 bc_mom_w_acc_ext=0.0,
6222 bc_dmom_w_acc_w_ext=0.0,
6223 bc_mass_adv_ext[nSpace],
6224 bc_dmass_adv_u_ext[nSpace],
6225 bc_dmass_adv_v_ext[nSpace],
6226 bc_dmass_adv_w_ext[nSpace],
6227 bc_mom_u_adv_ext[nSpace],
6228 bc_dmom_u_adv_u_ext[nSpace],
6229 bc_dmom_u_adv_v_ext[nSpace],
6230 bc_dmom_u_adv_w_ext[nSpace],
6231 bc_mom_v_adv_ext[nSpace],
6232 bc_dmom_v_adv_u_ext[nSpace],
6233 bc_dmom_v_adv_v_ext[nSpace],
6234 bc_dmom_v_adv_w_ext[nSpace],
6235 bc_mom_w_adv_ext[nSpace],
6236 bc_dmom_w_adv_u_ext[nSpace],
6237 bc_dmom_w_adv_v_ext[nSpace],
6238 bc_dmom_w_adv_w_ext[nSpace],
6239 bc_mom_uu_diff_ten_ext[nSpace],
6240 bc_mom_vv_diff_ten_ext[nSpace],
6241 bc_mom_ww_diff_ten_ext[nSpace],
6242 bc_mom_uv_diff_ten_ext[1],
6243 bc_mom_uw_diff_ten_ext[1],
6244 bc_mom_vu_diff_ten_ext[1],
6245 bc_mom_vw_diff_ten_ext[1],
6246 bc_mom_wu_diff_ten_ext[1],
6247 bc_mom_wv_diff_ten_ext[1],
6248 bc_mom_u_source_ext=0.0,
6249 bc_mom_v_source_ext=0.0,
6250 bc_mom_w_source_ext=0.0,
6251 bc_mom_u_ham_ext=0.0,
6252 bc_dmom_u_ham_grad_p_ext[nSpace],
6253 bc_dmom_u_ham_grad_u_ext[nSpace],
6254 bc_mom_v_ham_ext=0.0,
6255 bc_dmom_v_ham_grad_p_ext[nSpace],
6256 bc_dmom_v_ham_grad_v_ext[nSpace],
6257 bc_mom_w_ham_ext=0.0,
6258 bc_dmom_w_ham_grad_p_ext[nSpace],
6259 bc_dmom_w_ham_grad_w_ext[nSpace],
6260 fluxJacobian_p_p[nDOF_trial_element],
6261 fluxJacobian_p_u[nDOF_trial_element],
6262 fluxJacobian_p_v[nDOF_trial_element],
6263 fluxJacobian_p_w[nDOF_trial_element],
6264 fluxJacobian_u_p[nDOF_trial_element],
6265 fluxJacobian_u_u[nDOF_trial_element],
6266 fluxJacobian_u_v[nDOF_trial_element],
6267 fluxJacobian_u_w[nDOF_trial_element],
6268 fluxJacobian_v_p[nDOF_trial_element],
6269 fluxJacobian_v_u[nDOF_trial_element],
6270 fluxJacobian_v_v[nDOF_trial_element],
6271 fluxJacobian_v_w[nDOF_trial_element],
6272 fluxJacobian_w_p[nDOF_trial_element],
6273 fluxJacobian_w_u[nDOF_trial_element],
6274 fluxJacobian_w_v[nDOF_trial_element],
6275 fluxJacobian_w_w[nDOF_trial_element],
6276 jac_ext[nSpace*nSpace],
6278 jacInv_ext[nSpace*nSpace],
6279 boundaryJac[nSpace*(nSpace-1)],
6280 metricTensor[(nSpace-1)*(nSpace-1)],
6281 metricTensorDetSqrt,
6282 p_grad_trial_trace[nDOF_trial_element*nSpace],
6283 vel_grad_trial_trace[nDOF_trial_element*nSpace],
6285 p_test_dS[nDOF_test_element],
6286 vel_test_dS[nDOF_test_element],
6288 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
6289 vel_grad_test_dS[nDOF_trial_element*nSpace],
6293 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
6294 ck.calculateMapping_elementBoundary(eN,
6300 mesh_trial_trace_ref.data(),
6301 mesh_grad_trial_trace_ref.data(),
6302 boundaryJac_ref.data(),
6308 metricTensorDetSqrt,
6312 ck.calculateMappingVelocity_elementBoundary(eN,
6316 mesh_velocity_dof.data(),
6318 mesh_trial_trace_ref.data(),
6319 xt_ext,yt_ext,zt_ext,
6325 dS = metricTensorDetSqrt*dS_ref[kb];
6326 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
6327 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
6329 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
6330 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
6331 const double particle_eps = particle_epsFact * (useMetrics * h_phi + (1.0 - useMetrics) * elementDiameter[eN]);
6336 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
6339 p_ext = ebqe_p[ebNE_kb];
6340 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
6341 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
6344 for (
int I=0;I<nSpace;I++)
6345 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace+I];
6346 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
6347 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
6350 for (
int j=0;j<nDOF_trial_element;j++)
6353 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
6354 for (
int I=0;I<nSpace;I++)
6355 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
6360 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
6362 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
6363 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
6366 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
6370 double distance_to_omega_solid = 1e10;
6371 if (use_ball_as_particle == 1)
6373 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
6377 distance_to_omega_solid = ebq_global_phi_solid[ebN*nQuadraturePoints_elementBoundary+kb];
6379 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.),rhoSave, nuSave;
6388 elementDiameter[eN],
6389 smagorinskyConstant,
6390 turbulenceClosureModel,
6393 ebqe_vf_ext[ebNE_kb],
6394 ebqe_phi_ext[ebNE_kb],
6395 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
6396 distance_to_omega_solid,
6397 ebqe_kappa_phi_ext[ebNE_kb],
6409 ebqe_velocity_star[ebNE_kb_nSpace+0],
6410 ebqe_velocity_star[ebNE_kb_nSpace+1],
6411 ebqe_velocity_star[ebNE_kb_nSpace+1],
6435 mom_uu_diff_ten_ext,
6436 mom_vv_diff_ten_ext,
6437 mom_ww_diff_ten_ext,
6438 mom_uv_diff_ten_ext,
6439 mom_uw_diff_ten_ext,
6440 mom_vu_diff_ten_ext,
6441 mom_vw_diff_ten_ext,
6442 mom_wu_diff_ten_ext,
6443 mom_wv_diff_ten_ext,
6448 dmom_u_ham_grad_p_ext,
6449 dmom_u_ham_grad_u_ext,
6451 dmom_v_ham_grad_p_ext,
6452 dmom_v_ham_grad_v_ext,
6454 dmom_w_ham_grad_p_ext,
6455 dmom_w_ham_grad_w_ext,
6463 MATERIAL_PARAMETERS_AS_FUNCTION,
6464 ebqe_density_as_function[ebNE_kb],
6465 ebqe_dynamic_viscosity_as_function[ebNE_kb],
6468 use_ball_as_particle,
6471 ball_velocity.data(),
6472 ball_angular_velocity.data(),
6473 INT_BY_PARTS_PRESSURE);
6482 elementDiameter[eN],
6483 smagorinskyConstant,
6484 turbulenceClosureModel,
6487 bc_ebqe_vf_ext[ebNE_kb],
6488 bc_ebqe_phi_ext[ebNE_kb],
6489 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
6490 distance_to_omega_solid,
6491 ebqe_kappa_phi_ext[ebNE_kb],
6503 ebqe_velocity_star[ebNE_kb_nSpace+0],
6504 ebqe_velocity_star[ebNE_kb_nSpace+1],
6505 ebqe_velocity_star[ebNE_kb_nSpace+1],
6506 bc_eddy_viscosity_ext,
6508 bc_dmom_u_acc_u_ext,
6510 bc_dmom_v_acc_v_ext,
6512 bc_dmom_w_acc_w_ext,
6518 bc_dmom_u_adv_u_ext,
6519 bc_dmom_u_adv_v_ext,
6520 bc_dmom_u_adv_w_ext,
6522 bc_dmom_v_adv_u_ext,
6523 bc_dmom_v_adv_v_ext,
6524 bc_dmom_v_adv_w_ext,
6526 bc_dmom_w_adv_u_ext,
6527 bc_dmom_w_adv_v_ext,
6528 bc_dmom_w_adv_w_ext,
6529 bc_mom_uu_diff_ten_ext,
6530 bc_mom_vv_diff_ten_ext,
6531 bc_mom_ww_diff_ten_ext,
6532 bc_mom_uv_diff_ten_ext,
6533 bc_mom_uw_diff_ten_ext,
6534 bc_mom_vu_diff_ten_ext,
6535 bc_mom_vw_diff_ten_ext,
6536 bc_mom_wu_diff_ten_ext,
6537 bc_mom_wv_diff_ten_ext,
6538 bc_mom_u_source_ext,
6539 bc_mom_v_source_ext,
6540 bc_mom_w_source_ext,
6542 bc_dmom_u_ham_grad_p_ext,
6543 bc_dmom_u_ham_grad_u_ext,
6545 bc_dmom_v_ham_grad_p_ext,
6546 bc_dmom_v_ham_grad_v_ext,
6548 bc_dmom_w_ham_grad_p_ext,
6549 bc_dmom_w_ham_grad_w_ext,
6557 MATERIAL_PARAMETERS_AS_FUNCTION,
6558 ebqe_density_as_function[ebNE_kb],
6559 ebqe_dynamic_viscosity_as_function[ebNE_kb],
6562 use_ball_as_particle,
6565 ball_velocity.data(),
6566 ball_angular_velocity.data(),
6567 INT_BY_PARTS_PRESSURE);
6569 if (turbulenceClosureModel >= 3)
6571 const double turb_var_grad_0_dummy[2] = {0.,0.};
6572 const double c_mu = 0.09;
6581 ebqe_vf_ext[ebNE_kb],
6582 ebqe_phi_ext[ebNE_kb],
6585 ebqe_turb_var_0[ebNE_kb],
6586 ebqe_turb_var_1[ebNE_kb],
6587 turb_var_grad_0_dummy,
6589 mom_uu_diff_ten_ext,
6590 mom_vv_diff_ten_ext,
6591 mom_ww_diff_ten_ext,
6592 mom_uv_diff_ten_ext,
6593 mom_uw_diff_ten_ext,
6594 mom_vu_diff_ten_ext,
6595 mom_vw_diff_ten_ext,
6596 mom_wu_diff_ten_ext,
6597 mom_wv_diff_ten_ext,
6610 ebqe_vf_ext[ebNE_kb],
6611 ebqe_phi_ext[ebNE_kb],
6614 ebqe_turb_var_0[ebNE_kb],
6615 ebqe_turb_var_1[ebNE_kb],
6616 turb_var_grad_0_dummy,
6617 bc_eddy_viscosity_ext,
6618 bc_mom_uu_diff_ten_ext,
6619 bc_mom_vv_diff_ten_ext,
6620 bc_mom_ww_diff_ten_ext,
6621 bc_mom_uv_diff_ten_ext,
6622 bc_mom_uw_diff_ten_ext,
6623 bc_mom_vu_diff_ten_ext,
6624 bc_mom_vw_diff_ten_ext,
6625 bc_mom_wu_diff_ten_ext,
6626 bc_mom_wv_diff_ten_ext,
6627 bc_mom_u_source_ext,
6628 bc_mom_v_source_ext,
6629 bc_mom_w_source_ext);
6634 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
6635 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
6637 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
6638 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
6641 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
6642 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
6644 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
6645 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
6657 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
6658 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
6661 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
6662 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
6672 isDOFBoundary_u[ebNE_kb],
6673 isDOFBoundary_v[ebNE_kb],
6674 isDOFBoundary_w[ebNE_kb],
6675 isAdvectiveFluxBoundary_p[ebNE_kb],
6676 isAdvectiveFluxBoundary_u[ebNE_kb],
6677 isAdvectiveFluxBoundary_v[ebNE_kb],
6678 isAdvectiveFluxBoundary_w[ebNE_kb],
6679 dmom_u_ham_grad_p_ext[0],
6681 porosity_ext*dmom_u_acc_u_ext,
6690 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]),
6691 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
6692 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
6693 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
6720 dflux_mom_u_adv_p_ext,
6721 dflux_mom_u_adv_u_ext,
6722 dflux_mom_u_adv_v_ext,
6723 dflux_mom_u_adv_w_ext,
6724 dflux_mom_v_adv_p_ext,
6725 dflux_mom_v_adv_u_ext,
6726 dflux_mom_v_adv_v_ext,
6727 dflux_mom_v_adv_w_ext,
6728 dflux_mom_w_adv_p_ext,
6729 dflux_mom_w_adv_u_ext,
6730 dflux_mom_w_adv_v_ext,
6731 dflux_mom_w_adv_w_ext,
6732 &ebqe_velocity_star[ebNE_kb_nSpace]);
6736 ck.calculateGScale(G,normal,h_penalty);
6737 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
6738 for (
int j=0;j<nDOF_trial_element;j++)
6740 int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
6747 fluxJacobian_u_u[j] =
6748 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6750 ebqe_phi_ext[ebNE_kb],
6751 sdInfo_u_u_rowptr.data(),
6752 sdInfo_u_u_colind.data(),
6753 isDOFBoundary_u[ebNE_kb],
6754 isDiffusiveFluxBoundary_u[ebNE_kb],
6756 mom_uu_diff_ten_ext,
6757 vel_trial_trace_ref[ebN_local_kb_j],
6758 &vel_grad_trial_trace[j_nSpace],
6760 fluxJacobian_u_v[j]=
6761 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6763 ebqe_phi_ext[ebNE_kb],
6764 sdInfo_u_v_rowptr.data(),
6765 sdInfo_u_v_colind.data(),
6766 isDOFBoundary_v[ebNE_kb],
6767 isDiffusiveFluxBoundary_v[ebNE_kb],
6769 mom_uv_diff_ten_ext,
6770 vel_trial_trace_ref[ebN_local_kb_j],
6771 &vel_grad_trial_trace[j_nSpace],
6788 fluxJacobian_v_u[j]=
6789 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6791 ebqe_phi_ext[ebNE_kb],
6792 sdInfo_v_u_rowptr.data(),
6793 sdInfo_v_u_colind.data(),
6794 isDOFBoundary_u[ebNE_kb],
6795 isDiffusiveFluxBoundary_u[ebNE_kb],
6797 mom_vu_diff_ten_ext,
6798 vel_trial_trace_ref[ebN_local_kb_j],
6799 &vel_grad_trial_trace[j_nSpace],
6801 fluxJacobian_v_v[j]=
6802 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6804 ebqe_phi_ext[ebNE_kb],
6805 sdInfo_v_v_rowptr.data(),
6806 sdInfo_v_v_colind.data(),
6807 isDOFBoundary_v[ebNE_kb],
6808 isDiffusiveFluxBoundary_v[ebNE_kb],
6810 mom_vv_diff_ten_ext,
6811 vel_trial_trace_ref[ebN_local_kb_j],
6812 &vel_grad_trial_trace[j_nSpace],
6869 for (
int i=0;i<nDOF_test_element;i++)
6871 int eN_i = eN*nDOF_test_element+i;
6872 for (
int j=0;j<nDOF_trial_element;j++)
6874 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;
6882 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*vel_test_dS[i]+
6883 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6884 isDiffusiveFluxBoundary_u[ebNE_kb],
6886 vel_trial_trace_ref[ebN_local_kb_j],
6888 sdInfo_u_u_rowptr.data(),
6889 sdInfo_u_u_colind.data(),
6890 mom_uu_diff_ten_ext,
6891 &vel_grad_test_dS[i*nSpace]);
6892 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] += fluxJacobian_u_v[j]*vel_test_dS[i]+
6893 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6894 isDiffusiveFluxBoundary_u[ebNE_kb],
6896 vel_trial_trace_ref[ebN_local_kb_j],
6898 sdInfo_u_v_rowptr.data(),
6899 sdInfo_u_v_colind.data(),
6900 mom_uv_diff_ten_ext,
6901 &vel_grad_test_dS[i*nSpace]);
6914 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] += fluxJacobian_v_u[j]*vel_test_dS[i]+
6915 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6916 isDiffusiveFluxBoundary_v[ebNE_kb],
6918 vel_trial_trace_ref[ebN_local_kb_j],
6920 sdInfo_v_u_rowptr.data(),
6921 sdInfo_v_u_colind.data(),
6922 mom_vu_diff_ten_ext,
6923 &vel_grad_test_dS[i*nSpace]);
6924 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] += fluxJacobian_v_v[j]*vel_test_dS[i]+
6925 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6926 isDiffusiveFluxBoundary_v[ebNE_kb],
6928 vel_trial_trace_ref[ebN_local_kb_j],
6930 sdInfo_v_v_rowptr.data(),
6931 sdInfo_v_v_colind.data(),
6932 mom_vv_diff_ten_ext,
6933 &vel_grad_test_dS[i*nSpace]);