25 #include "libmesh/coupling_matrix.h" 26 #include "libmesh/dof_map.h" 27 #include "libmesh/elem.h" 28 #include "libmesh/equation_systems.h" 29 #include "libmesh/fe_interface.h" 30 #include "libmesh/node.h" 31 #include "libmesh/quadrature_gauss.h" 32 #include "libmesh/sparse_matrix.h" 33 #include "libmesh/tensor_value.h" 34 #include "libmesh/vector_value.h" 35 #include "libmesh/fe.h" 37 template <
typename P,
typename C>
48 template <
typename P,
typename C>
57 ?
mesh.getCoordSystem(sub_id) ==
mesh.getCoordSystem(neighbor_sub_id)
59 "Coordinate systems must be the same between element and neighbor");
60 const auto coord_type =
mesh.getCoordSystem(sub_id);
64 if (
mesh.usingGeneralAxisymmetricCoordAxes())
66 const auto & axis =
mesh.getGeneralAxisymmetricCoordAxis(sub_id);
71 point, factor, coord_type,
mesh.getAxisymmetricRadialCoord());
79 _subproblem(_sys.subproblem()),
81 _nonlocal_cm(_subproblem.nonlocalCouplingMatrix(_sys.number())),
82 _computing_residual(_subproblem.currentlyComputingResidual()),
83 _computing_jacobian(_subproblem.currentlyComputingJacobian()),
84 _computing_residual_and_jacobian(_subproblem.currentlyComputingResidualAndJacobian()),
85 _dof_map(_sys.dofMap()),
88 _mesh_dimension(_mesh.dimension()),
90 _user_added_fe_of_helper_type(false),
91 _user_added_fe_face_of_helper_type(false),
92 _user_added_fe_face_neighbor_of_helper_type(false),
93 _user_added_fe_neighbor_of_helper_type(false),
94 _user_added_fe_lower_of_helper_type(false),
95 _building_helpers(false),
96 _current_qrule(nullptr),
97 _current_qrule_volume(nullptr),
98 _current_qrule_arbitrary(nullptr),
100 _current_qrule_face(nullptr),
101 _current_qface_arbitrary(nullptr),
102 _current_qrule_neighbor(nullptr),
103 _need_JxW_neighbor(false),
105 _custom_mortar_qrule(false),
106 _current_qrule_lower(nullptr),
108 _current_elem(nullptr),
109 _current_elem_volume(0),
111 _current_side_elem(nullptr),
112 _current_side_volume(0),
113 _current_neighbor_elem(nullptr),
114 _current_neighbor_side(0),
115 _current_neighbor_side_elem(nullptr),
116 _need_neighbor_elem_volume(false),
117 _current_neighbor_volume(0),
118 _current_node(nullptr),
119 _current_neighbor_node(nullptr),
120 _current_elem_volume_computed(false),
121 _current_side_volume_computed(false),
123 _current_lower_d_elem(nullptr),
124 _current_neighbor_lower_d_elem(nullptr),
125 _need_lower_d_elem_volume(false),
126 _need_neighbor_lower_d_elem_volume(false),
130 _cached_residual_values(2),
131 _cached_residual_rows(2),
132 _max_cached_residuals(0),
133 _max_cached_jacobians(0),
135 _block_diagonal_matrix(false),
136 _calculate_xyz(false),
137 _calculate_face_xyz(false),
138 _calculate_curvatures(false),
139 _calculate_ad_coord(false),
140 _have_p_refinement(false)
174 _fe_msm->add_p_level_in_reinit(
false);
187 for (
auto & it :
_fe[
dim])
276 _fe[
dim][type] = FEGenericBase<Real>::build(
dim, type).release();
278 _fe[
dim][type]->get_phi();
279 _fe[
dim][type]->get_dphi();
283 _fe[
dim][type]->get_xyz();
285 _fe[
dim][type]->get_d2phi();
302 _fe_face[
dim][type] = FEGenericBase<Real>::build(
dim, type).release();
370 _fe_lower[
dim][type] = FEGenericBase<Real>::build(
dim, type).release();
391 _fe_lower[
dim][type] = FEGenericBase<Real>::build(
dim, type).release();
411 if (ending_dim <
dim)
413 for (;
dim <= ending_dim;
dim++)
436 if (ending_dim <
dim)
438 for (;
dim <= ending_dim;
dim++)
457 unsigned int min_dim;
468 _vector_fe[
dim][type] = FEGenericBase<VectorValue<Real>>::build(
dim, type).release();
490 unsigned int min_dim;
519 unsigned int min_dim;
548 unsigned int min_dim;
560 FEGenericBase<VectorValue<Real>>::build(
dim, type).release();
575 mooseAssert(qdefault.size() > 0,
"default quadrature must be initialized before order bumps");
579 if (qvec.size() != ndims || !qvec[0].vol)
581 qdefault[0].arbitrary_vol->get_order(),
583 qdefault[0].face->get_order(),
585 else if (qvec[0].vol->get_order() < volume_order)
587 qvec[0].arbitrary_vol->get_order(),
589 qvec[0].face->get_order(),
598 mooseAssert(qdefault.size() > 0,
"default quadrature must be initialized before order bumps");
602 if (qvec.size() != ndims || !qvec[0].vol)
603 createQRules(qdefault[0].vol->type(), order, order, order, block);
604 else if (qvec[0].vol->get_order() < order || qvec[0].face->get_order() < order)
606 std::max(order, qvec[0].arbitrary_vol->get_order()),
607 std::max(order, qvec[0].vol->get_order()),
608 std::max(order, qvec[0].face->get_order()),
619 bool allow_negative_qweights)
623 if (qvec.size() != ndims)
626 for (
unsigned int i = 0; i < qvec.size(); i++)
629 auto & q = qvec[
dim];
630 q.vol = QBase::build(type,
dim, volume_order);
631 q.vol->allow_rules_with_negative_weights = allow_negative_qweights;
632 q.face = QBase::build(type,
dim - 1, face_order);
633 q.face->allow_rules_with_negative_weights = allow_negative_qweights;
635 q.fv_face->allow_rules_with_negative_weights = allow_negative_qweights;
636 q.neighbor = std::make_unique<ArbitraryQuadrature>(
dim - 1, face_order);
637 q.neighbor->allow_rules_with_negative_weights = allow_negative_qweights;
638 q.arbitrary_vol = std::make_unique<ArbitraryQuadrature>(
dim, order);
639 q.arbitrary_vol->allow_rules_with_negative_weights = allow_negative_qweights;
640 q.arbitrary_face = std::make_unique<ArbitraryQuadrature>(
dim - 1, face_order);
641 q.arbitrary_face->allow_rules_with_negative_weights = allow_negative_qweights;
647 _qrule_msm->allow_rules_with_negative_weights = allow_negative_qweights;
658 for (
auto & it :
_fe[
dim])
659 it.second->attach_quadrature_rule(qrule);
661 it.second->attach_quadrature_rule(qrule);
676 it.second->attach_quadrature_rule(qrule);
678 it.second->attach_quadrature_rule(qrule);
695 it.second->attach_quadrature_rule(qrule);
697 it.second->attach_quadrature_rule(qrule);
711 it.second->attach_quadrature_rule(qrule);
713 it.second->attach_quadrature_rule(qrule);
717 "We should not be indexing out of bounds");
750 " does not match previously specified quadrature_order: ",
752 ". Quadrature_order (when specified) must match for all mortar constraints.");
759 unsigned int dim =
elem->dim();
761 for (
const auto & it :
_fe[
dim])
764 const FEType & fe_type = it.first;
772 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe.get_phi()));
774 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe.get_dphi()));
777 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe.get_d2phi()));
782 const FEType & fe_type = it.first;
790 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe.get_phi()));
792 const_cast<std::vector<std::vector<TensorValue<Real>
>> &>(fe.get_dphi()));
795 const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe.get_d2phi()));
800 fesd.
_div_phi.shallowCopy(const_cast<std::vector<std::vector<Real>> &>(fe.get_div_phi()));
821 for (
unsigned int qp = 0; qp != n_qp; qp++)
825 for (
unsigned qp = 0; qp < n_qp; ++qp)
832 for (
const auto & it :
_fe[
dim])
835 auto fe_type = it.first;
836 auto num_shapes = fe.n_shape_functions();
839 grad_phi.resize(num_shapes);
840 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
841 grad_phi[i].resize(n_qp);
847 const auto & regular_grad_phi =
_fe_shape_data[fe_type]->_grad_phi;
848 for (
unsigned qp = 0; qp < n_qp; ++qp)
849 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
850 grad_phi[i][qp] = regular_grad_phi[i][qp];
856 auto fe_type = it.first;
857 auto num_shapes = fe.n_shape_functions();
860 grad_phi.resize(num_shapes);
861 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
862 grad_phi[i].resize(n_qp);
869 for (
unsigned qp = 0; qp < n_qp; ++qp)
870 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
871 grad_phi[i][qp] = regular_grad_phi[i][qp];
881 if (
_xfem !=
nullptr)
885 template <
typename OutputType>
890 FEGenericBase<OutputType> * fe)
904 const auto & dphidxi = fe->get_dphidxi();
905 const auto & dphideta = fe->get_dphideta();
906 const auto & dphidzeta = fe->get_dphidzeta();
907 auto num_shapes = grad_phi.size();
913 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
914 for (
unsigned qp = 0; qp < n_qp; ++qp)
921 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
922 for (
unsigned qp = 0; qp < n_qp; ++qp)
924 grad_phi[i][qp].slice(0) = dphidxi[i][qp] *
_ad_dxidx_map[qp];
925 grad_phi[i][qp].slice(1) = dphidxi[i][qp] *
_ad_dxidy_map[qp];
926 grad_phi[i][qp].slice(2) = dphidxi[i][qp] *
_ad_dxidz_map[qp];
933 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
934 for (
unsigned qp = 0; qp < n_qp; ++qp)
936 grad_phi[i][qp].slice(0) =
938 grad_phi[i][qp].slice(1) =
940 grad_phi[i][qp].slice(2) =
948 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
949 for (
unsigned qp = 0; qp < n_qp; ++qp)
951 grad_phi[i][qp].slice(0) = dphidxi[i][qp] *
_ad_dxidx_map[qp] +
954 grad_phi[i][qp].slice(1) = dphidxi[i][qp] *
_ad_dxidy_map[qp] +
957 grad_phi[i][qp].slice(2) = dphidxi[i][qp] *
_ad_dxidz_map[qp] +
998 const std::vector<Real> & qw,
1040 const auto & elem_nodes =
elem->get_nodes();
1041 auto num_shapes = fe->n_shape_functions();
1042 const auto & phi_map = fe->get_fe_map().get_phi_map();
1043 const auto & dphidxi_map = fe->get_fe_map().get_dphidxi_map();
1044 const auto & dphideta_map = fe->get_fe_map().get_dphideta_map();
1045 const auto & dphidzeta_map = fe->get_fe_map().get_dphidzeta_map();
1047 const bool do_derivatives =
1068 for (std::size_t i = 0; i < num_shapes; i++)
1071 const Node &
node = *elem_nodes[i];
1075 if (
node.n_dofs(sys_num, disp_num))
1077 elem_point(direction).derivatives(),
node.dof_number(sys_num, disp_num, 0), 1.);
1089 static bool failing =
false;
1094 libmesh_error_msg(
"ERROR: negative Jacobian " <<
_ad_jac[p].
value() <<
" at point index " 1095 << p <<
" in element " <<
elem->id());
1118 for (std::size_t i = 0; i < num_shapes; i++)
1121 const Node &
node = *elem_nodes[i];
1126 elem_point(direction).derivatives(),
node.dof_number(sys_num, disp_num, 0), 1.);
1139 const auto g11 = (dx_dxi * dx_dxi + dy_dxi * dy_dxi + dz_dxi * dz_dxi);
1141 const auto g12 = (dx_dxi * dx_deta + dy_dxi * dy_deta + dz_dxi * dz_deta);
1143 const auto g21 = g12;
1145 const auto g22 = (dx_deta * dx_deta + dy_deta * dy_deta + dz_deta * dz_deta);
1147 auto det = (g11 * g22 - g12 * g21);
1149 if (det.value() <= -TOLERANCE * TOLERANCE)
1151 static bool failing =
false;
1156 libmesh_error_msg(
"ERROR: negative Jacobian " << det <<
" at point index " << p
1157 <<
" in element " <<
elem->id());
1162 else if (det.value() <= 0.)
1163 det.value() = TOLERANCE * TOLERANCE;
1165 const auto inv_det = 1. / det;
1170 const auto g11inv = g22 * inv_det;
1171 const auto g12inv = -g12 * inv_det;
1172 const auto g21inv = -g21 * inv_det;
1173 const auto g22inv = g11 * inv_det;
1194 for (std::size_t i = 0; i < num_shapes; i++)
1197 const Node &
node = *elem_nodes[i];
1202 elem_point(direction).derivatives(),
node.dof_number(sys_num, disp_num, 0), 1.);
1218 _ad_jac[p] = (dx_dxi * (dy_deta * dz_dzeta - dz_deta * dy_dzeta) +
1219 dy_dxi * (dz_deta * dx_dzeta - dx_deta * dz_dzeta) +
1220 dz_dxi * (dx_deta * dy_dzeta - dy_deta * dx_dzeta));
1224 static bool failing =
false;
1229 libmesh_error_msg(
"ERROR: negative Jacobian " <<
_ad_jac[p].
value() <<
" at point index " 1230 << p <<
" in element " <<
elem->id());
1238 const auto inv_jac = 1. /
_ad_jac[p];
1240 _ad_dxidx_map[p] = (dy_deta * dz_dzeta - dz_deta * dy_dzeta) * inv_jac;
1241 _ad_dxidy_map[p] = (dz_deta * dx_dzeta - dx_deta * dz_dzeta) * inv_jac;
1242 _ad_dxidz_map[p] = (dx_deta * dy_dzeta - dy_deta * dx_dzeta) * inv_jac;
1244 _ad_detadx_map[p] = (dz_dxi * dy_dzeta - dy_dxi * dz_dzeta) * inv_jac;
1245 _ad_detady_map[p] = (dx_dxi * dz_dzeta - dz_dxi * dx_dzeta) * inv_jac;
1246 _ad_detadz_map[p] = (dy_dxi * dx_dzeta - dx_dxi * dy_dzeta) * inv_jac;
1256 libmesh_error_msg(
"Invalid dim = " <<
dim);
1263 unsigned int dim =
elem->dim();
1267 FEBase & fe_face = *it.second;
1268 const FEType & fe_type = it.first;
1273 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_face.get_phi()));
1275 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe_face.get_dphi()));
1278 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
1283 const FEType & fe_type = it.first;
1291 fesd.
_phi.shallowCopy(
1292 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe_face.get_phi()));
1294 const_cast<std::vector<std::vector<TensorValue<Real>
>> &>(fe_face.get_dphi()));
1297 const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
1300 const_cast<std::vector<std::vector<
VectorValue<Real>>> &>(fe_face.get_curl_phi()));
1303 const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
1331 if (
_xfem !=
nullptr)
1351 const auto n_qp = qw.size();
1355 std::vector<std::vector<Real>>
const * d2psidxi2_map =
nullptr;
1356 std::vector<std::vector<Real>>
const * d2psidxideta_map =
nullptr;
1357 std::vector<std::vector<Real>>
const * d2psideta2_map =
nullptr;
1375 if (side_elem.node_id(0) ==
elem.node_id(0))
1380 VectorValue<DualReal> side_point;
1383 const Node &
node = side_elem.node_ref(0);
1389 side_point(direction).derivatives(),
node.dof_number(sys_num, disp_num, 0), 1.);
1392 for (
unsigned int p = 0; p < n_qp; p++)
1413 for (
unsigned int p = 0; p < n_qp; p++)
1422 const auto n_mapping_shape_functions =
1423 FE<2, LAGRANGE>::n_shape_functions(side_elem.type(), side_elem.default_order());
1425 for (
unsigned int i = 0; i < n_mapping_shape_functions; i++)
1427 const Node &
node = side_elem.node_ref(i);
1428 VectorValue<DualReal> side_point =
node;
1433 side_point(direction).derivatives(),
node.dof_number(sys_num, disp_num, 0), 1.);
1435 for (
unsigned int p = 0; p < n_qp; p++)
1445 for (
unsigned int p = 0; p < n_qp; p++)
1455 libmesh_assert_not_equal_to(denominator, 0);
1474 for (
unsigned int p = 0; p < n_qp; p++)
1488 const unsigned int n_mapping_shape_functions =
1489 FE<3, LAGRANGE>::n_shape_functions(side_elem.type(), side_elem.default_order());
1491 for (
unsigned int i = 0; i < n_mapping_shape_functions; i++)
1493 const Node &
node = side_elem.node_ref(i);
1494 VectorValue<DualReal> side_point =
node;
1499 side_point(direction).derivatives(),
node.dof_number(sys_num, disp_num, 0), 1.);
1501 for (
unsigned int p = 0; p < n_qp; p++)
1516 for (
unsigned int p = 0; p < n_qp; p++)
1524 const auto g11 = (dxdxi * dxdxi + dydxi * dydxi + dzdxi * dzdxi);
1526 const auto g12 = (dxdxi * dxdeta + dydxi * dydeta + dzdxi * dzdeta);
1528 const auto g21 = g12;
1530 const auto g22 = (dxdeta * dxdeta + dydeta * dydeta + dzdeta * dzdeta);
1532 const auto the_jac =
std::sqrt(g11 * g22 - g12 * g21);
1545 const auto numerator = E *
N - 2. * F * M + G * L;
1546 const auto denominator = E * G - F * F;
1547 libmesh_assert_not_equal_to(denominator, 0.);
1563 unsigned int neighbor_dim =
neighbor->dim();
1568 FEBase & fe_face_neighbor = *it.second;
1569 FEType fe_type = it.first;
1572 fe_face_neighbor.reinit(
neighbor, &reference_points);
1576 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_face_neighbor.get_phi()));
1578 const_cast<std::vector<std::vector<RealGradient>
> &>(fe_face_neighbor.get_dphi()));
1581 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
1586 const FEType & fe_type = it.first;
1592 fe_face_neighbor.reinit(
neighbor, &reference_points);
1594 fesd.
_phi.shallowCopy(
1595 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe_face_neighbor.get_phi()));
1597 const_cast<std::vector<std::vector<TensorValue<Real>
>> &>(fe_face_neighbor.get_dphi()));
1599 fesd.
_second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
1600 fe_face_neighbor.get_d2phi()));
1603 fe_face_neighbor.get_curl_phi()));
1606 const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
1611 "We should be in bounds here");
1622 unsigned int neighbor_dim =
neighbor->dim();
1627 FEBase & fe_neighbor = *it.second;
1628 FEType fe_type = it.first;
1631 fe_neighbor.reinit(
neighbor, &reference_points);
1635 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_neighbor.get_phi()));
1637 const_cast<std::vector<std::vector<RealGradient>
> &>(fe_neighbor.get_dphi()));
1640 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_neighbor.get_d2phi()));
1645 const FEType & fe_type = it.first;
1651 fe_neighbor.reinit(
neighbor, &reference_points);
1653 fesd.
_phi.shallowCopy(
1654 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe_neighbor.get_phi()));
1656 const_cast<std::vector<std::vector<TensorValue<Real>
>> &>(fe_neighbor.get_dphi()));
1659 const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_neighbor.get_d2phi()));
1662 const_cast<std::vector<std::vector<
VectorValue<Real>>> &>(fe_neighbor.get_curl_phi()));
1665 const_cast<std::vector<std::vector<Real>> &>(fe_neighbor.get_div_phi()));
1677 unsigned int neighbor_dim =
neighbor->dim();
1679 "Neighbor subdomain ID has not been correctly set");
1683 neighbor_rule->setPoints(reference_points);
1688 "current neighbor subdomain has been set incorrectly");
1697 fe.attach_quadrature_rule(qrule);
1700 const std::vector<Real> &
JxW = fe.get_JxW();
1702 q_points.
shallowCopy(
const_cast<std::vector<Point> &
>(fe.get_xyz()));
1707 for (
unsigned int qp = 0; qp < qrule->n_points(); qp++)
1717 template <
typename Po
ints,
typename Coords>
1720 const Points & q_points,
1725 mooseAssert(qrule,
"The quadrature rule is null in Assembly::setCoordinateTransformation");
1726 auto n_points = qrule->n_points();
1727 mooseAssert(n_points == q_points.size(),
1728 "The number of points in the quadrature rule doesn't match the number of passed-in " 1729 "points in Assembly::setCoordinateTransformation");
1736 coord.resize(n_points);
1737 for (
unsigned int qp = 0; qp < n_points; qp++)
1785 "current subdomain has been set incorrectly");
1788 FEInterface::inverse_map(
elem->dim(),
1803 unsigned int elem_dimension =
elem->dim();
1816 "current subdomain has been set incorrectly");
1831 "current subdomain has been set incorrectly");
1857 "current subdomain has been set incorrectly");
1881 "Our finite volume quadrature rule should always yield a single point");
1889 const auto & ref_point = ref_points[0];
1896 "current neighbor subdomain has been set incorrectly");
1902 neighbor_rule->setPoints(ref_points);
1912 return qruleFaceHelper<QBase>(
elem,
side, [](
QRules & q) {
return q.face.get(); });
1918 return qruleFaceHelper<ArbitraryQuadrature>(
1925 const auto elem_dimension =
elem->dim();
1938 "current subdomain has been set incorrectly");
1952 Assembly::reinit(
const Elem * elem,
unsigned int side,
const std::vector<Point> & reference_points)
1957 "current subdomain has been set incorrectly");
1989 const Elem * neighbor,
1990 unsigned int neighbor_side,
1991 const std::vector<Point> * neighbor_reference_points)
1997 unsigned int neighbor_dim =
neighbor->dim();
1999 if (neighbor_reference_points)
2002 FEInterface::inverse_map(neighbor_dim,
2016 unsigned int elem_side,
2018 const std::vector<Point> *
const pts,
2019 const std::vector<Real> *
const weights)
2023 unsigned int elem_dim =
elem->dim();
2029 face_rule->setPoints(*pts);
2040 for (
const auto & it :
_fe_face[elem_dim])
2042 FEBase & fe_face = *it.second;
2043 FEType fe_type = it.first;
2046 fe_face.reinit(
elem, elem_side, tolerance, pts, weights);
2050 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_face.get_phi()));
2052 const_cast<std::vector<std::vector<RealGradient>
> &>(fe_face.get_dphi()));
2055 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face.get_d2phi()));
2060 const FEType & fe_type = it.first;
2066 fe_face.reinit(
elem, elem_side, tolerance, pts, weights);
2068 fesd.
_phi.shallowCopy(
2069 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe_face.get_phi()));
2071 const_cast<std::vector<std::vector<TensorValue<Real>
>> &>(fe_face.get_dphi()));
2074 const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(fe_face.get_d2phi()));
2077 const_cast<std::vector<std::vector<
VectorValue<Real>>> &>(fe_face.get_curl_phi()));
2080 const_cast<std::vector<std::vector<Real>> &>(fe_face.get_div_phi()));
2126 const std::vector<Real> dummy_qw(n_qp, 1.);
2128 for (
unsigned int qp = 0; qp != n_qp; qp++)
2132 for (
unsigned qp = 0; qp < n_qp; ++qp)
2144 FEBase & fe = *it.second;
2145 auto fe_type = it.first;
2146 auto num_shapes = fe.n_shape_functions();
2149 grad_phi.resize(num_shapes);
2150 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2151 grad_phi[i].resize(n_qp);
2158 for (
unsigned qp = 0; qp < n_qp; ++qp)
2159 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2160 grad_phi[i][qp] = regular_grad_phi[i][qp];
2165 auto fe_type = it.first;
2166 auto num_shapes = fe.n_shape_functions();
2169 grad_phi.resize(num_shapes);
2170 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2171 grad_phi[i].resize(n_qp);
2178 for (
unsigned qp = 0; qp < n_qp; ++qp)
2179 for (decltype(num_shapes) i = 0; i < num_shapes; ++i)
2180 grad_phi[i][qp] = regular_grad_phi[i][qp];
2187 unsigned int neighbor_side,
2189 const std::vector<Point> *
const pts,
2190 const std::vector<Real> *
const weights)
2194 unsigned int neighbor_dim =
neighbor->dim();
2210 FEBase & fe_face_neighbor = *it.second;
2211 FEType fe_type = it.first;
2214 fe_face_neighbor.reinit(
neighbor, neighbor_side, tolerance, pts, weights);
2218 fesd.
_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_face_neighbor.get_phi()));
2220 const_cast<std::vector<std::vector<RealGradient>
> &>(fe_face_neighbor.get_dphi()));
2223 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_face_neighbor.get_d2phi()));
2228 const FEType & fe_type = it.first;
2234 fe_face_neighbor.reinit(
neighbor, neighbor_side, tolerance, pts, weights);
2236 fesd.
_phi.shallowCopy(
2237 const_cast<std::vector<std::vector<VectorValue<Real>
>> &>(fe_face_neighbor.get_phi()));
2239 const_cast<std::vector<std::vector<TensorValue<Real>
>> &>(fe_face_neighbor.get_dphi()));
2241 fesd.
_second_phi.shallowCopy(const_cast<std::vector<std::vector<TypeNTensor<3, Real>>> &>(
2242 fe_face_neighbor.get_d2phi()));
2245 fe_face_neighbor.get_curl_phi()));
2248 const_cast<std::vector<std::vector<Real>> &>(fe_face_neighbor.get_div_phi()));
2253 "We should be in bounds here");
2255 neighbor, neighbor_side, tolerance, pts, weights);
2265 const std::vector<Point> & pts,
2266 const std::vector<Real> & JxW)
2268 const unsigned int elem_dim =
elem->dim();
2270 "Dual shape functions should only be computed on lower dimensional face elements");
2272 for (
const auto & it :
_fe_lower[elem_dim])
2274 FEBase & fe_lower = *it.second;
2276 fe_lower.set_calculate_default_dual_coeff(
false);
2277 fe_lower.reinit_dual_shape_coeffs(
elem, pts,
JxW);
2283 const std::vector<Point> *
const pts,
2284 const std::vector<Real> *
const weights)
2288 const unsigned int elem_dim =
elem->dim();
2290 "The lower dimensional element should truly be a lower dimensional element");
2308 for (
const auto & it :
_fe_lower[elem_dim])
2310 FEBase & fe_lower = *it.second;
2311 FEType fe_type = it.first;
2313 fe_lower.reinit(
elem);
2317 fesd->_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_lower.get_phi()));
2318 fesd->_grad_phi.shallowCopy(
2319 const_cast<std::vector<std::vector<RealGradient>
> &>(fe_lower.get_dphi()));
2321 fesd->_second_phi.shallowCopy(
2322 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_d2phi()));
2328 fesd->_phi.shallowCopy(
const_cast<std::vector<std::vector<Real>
> &>(fe_lower.get_dual_phi()));
2329 fesd->_grad_phi.shallowCopy(
2330 const_cast<std::vector<std::vector<RealGradient>
> &>(fe_lower.get_dual_dphi()));
2332 fesd->_second_phi.shallowCopy(
2333 const_cast<std::vector<std::vector<TensorValue<Real>>> &>(fe_lower.get_dual_d2phi()));
2345 if (pts && !weights)
2363 const auto & physical_q_points = helper_fe.get_xyz();
2364 const auto &
JxW = helper_fe.get_JxW();
2378 "You should be calling reinitNeighborLowerDElem on a lower dimensional element");
2399 "You should be calling reinitMortarElem on a lower dimensional element");
2411 unsigned int neighbor_side,
2412 const std::vector<Point> & physical_points)
2414 unsigned int neighbor_dim =
neighbor->dim();
2415 FEInterface::inverse_map(
2421 physical_points.size() == 1,
2422 "If reinitializing with more than one point, then I am dubious of your use case. Perhaps " 2423 "you are performing a DG type method and you are reinitializing using points from the " 2424 "element face. In such a case your neighbor JxW must have its index order 'match' the " 2425 "element JxW index order, e.g. imagining a vertical 1D face with two quadrature points, " 2427 "index 0 for elem JxW corresponds to the 'top' quadrature point, then index 0 for " 2429 "JxW must also correspond to the 'top' quadrature point. And libMesh/MOOSE has no way to " 2430 "guarantee that with multiple quadrature points.");
2448 const std::vector<Point> & physical_points)
2450 unsigned int neighbor_dim =
neighbor->dim();
2451 FEInterface::inverse_map(
2475 for (
auto & ivar : vars)
2477 auto i = ivar->number();
2482 for (
unsigned int k = 0; k < ivar->count(); ++k)
2484 unsigned int iv = i + k;
2485 for (
const auto & j : ConstCouplingRow(iv, *
_cm))
2496 auto pair = std::make_pair(ivar, &jvar);
2497 auto c = ivar_start;
2499 bool has_pair =
false;
2512 if (i != jvar.number())
2523 for (
auto & ivar : scalar_vars)
2525 auto i = ivar->number();
2529 for (
const auto & j : ConstCouplingRow(i, *
_cm))
2547 _sub_Re.resize(num_vector_tags);
2548 _sub_Rn.resize(num_vector_tags);
2549 _sub_Rl.resize(num_vector_tags);
2583 for (MooseIndex(num_matrix_tags) tag = 0; tag < num_matrix_tags; tag++)
2649 for (
auto & ivar : vars)
2651 auto i = ivar->number();
2653 for (
unsigned int k = 0; k < ivar->count(); ++k)
2655 unsigned int iv = i + k;
2656 for (
const auto & j : ConstCouplingRow(iv,
_nonlocal_cm))
2660 auto pair = std::make_pair(ivar, &jvar);
2661 auto c = ivar_start;
2663 bool has_pair =
false;
2685 unsigned int vi = ivar.
number();
2686 unsigned int vj = jvar.
number();
2703 for (
const auto & var : vars)
2705 tag_Re[var->number()].resize(var->dofIndices().size() * var->count());
2723 unsigned int vi = ivar.
number();
2724 unsigned int vj = jvar.
number();
2747 unsigned int vi = ivar.
number();
2748 unsigned int vj = jvar.
number();
2775 unsigned int vi = ivar.
number();
2776 unsigned int vj = jvar.
number();
2802 unsigned int vi = ivar.
number();
2803 unsigned int vj = jvar.
number();
2828 for (
const auto & var : vars)
2830 tag_Rn[var->number()].resize(var->dofIndicesNeighbor().size() * var->count());
2841 unsigned int vi = ivar.
number();
2842 unsigned int vj = jvar.
number();
2884 for (
const auto & var : vars)
2886 tag_Rl[var->number()].resize(var->dofIndicesLower().size() * var->count());
2892 const std::vector<dof_id_type> & dof_indices)
2896 const unsigned int ivn = iv.
number();
2897 const unsigned int jvn = jv.number();
2898 const unsigned int icount = iv.count();
2899 unsigned int jcount = jv.count();
2906 .resize(dof_indices.size() * icount, dof_indices.size() * jcount);
2911 tag_Re[ivn].resize(dof_indices.size() * icount);
2917 const std::vector<dof_id_type> & idof_indices,
2918 const std::vector<dof_id_type> & jdof_indices)
2922 const unsigned int ivn = iv.
number();
2923 const unsigned int jvn = jv.number();
2924 const unsigned int icount = iv.count();
2925 unsigned int jcount = jv.count();
2934 .resize(idof_indices.size() * icount, jdof_indices.size() * jcount);
2944 for (
const auto & ivar : vars)
2946 auto idofs = ivar->dofIndices().size();
2949 tag_Re[ivar->number()].resize(idofs);
2951 for (
const auto & jvar : vars)
2953 auto jdofs = jvar->dofIndices().size();
2970 for (
const auto & ivar : scalar_vars)
2972 auto idofs = ivar->dofIndices().size();
2974 for (
const auto & jvar : vars)
2976 auto jdofs = jvar->dofIndices().size() * jvar->count();
2989 template <
typename T>
2993 phi(v).shallowCopy(v.
phi());
3017 if (v.computingCurl())
3018 curlPhi(v).shallowCopy(v.curlPhi());
3019 if (v.computingDiv())
3020 divPhi(v).shallowCopy(v.divPhi());
3023 mooseError(
"Unsupported variable field type!");
3026 template <
typename T>
3054 if (v.computingCurl())
3056 if (v.computingDiv())
3060 mooseError(
"Unsupported variable field type!");
3063 template <
typename T>
3104 mooseError(
"Unsupported variable field type!");
3210 std::vector<dof_id_type> & dof_indices,
3211 const std::vector<Real> & scaling_factor,
3217 auto ndof = dof_indices.size();
3218 auto ntdof = res_block.
size();
3219 auto count = ntdof / ndof;
3220 mooseAssert(count == scaling_factor.size(),
"Inconsistent of number of components");
3221 mooseAssert(count * ndof == ntdof,
"Inconsistent of number of components");
3225 dof_indices.resize(ntdof);
3227 for (MooseIndex(count) j = 0; j < count; ++j)
3228 for (MooseIndex(ndof) i = 0; i < ndof; ++i)
3230 dof_indices[p] = dof_indices[i] + (is_nodal ? j : j * ndof);
3231 res_block(p) *= scaling_factor[j];
3237 if (scaling_factor[0] != 1.0)
3238 res_block *= scaling_factor[0];
3241 _dof_map.constrain_element_vector(res_block, dof_indices,
false);
3247 const std::vector<dof_id_type> & dof_indices,
3248 const std::vector<Real> & scaling_factor,
3251 if (dof_indices.size() > 0 && res_block.
size())
3262 std::vector<dof_id_type> & cached_residual_rows,
3264 const std::vector<dof_id_type> & dof_indices,
3265 const std::vector<Real> & scaling_factor,
3268 if (dof_indices.size() > 0 && res_block.
size())
3276 cached_residual_values.push_back(
_tmp_Re(i));
3287 const std::vector<dof_id_type> & dof_indices,
3288 const std::vector<Real> & scaling_factor,
3291 if (dof_indices.size() > 0)
3293 std::vector<dof_id_type> di(dof_indices);
3304 "Non-residual tag in Assembly::addResidual");
3309 for (
const auto & var : vars)
3311 tag_Re[var->number()],
3313 var->arrayScalingFactor(),
3320 for (
const auto & vector_tag : vector_tags)
3329 "Non-residual tag in Assembly::addResidualNeighbor");
3334 for (
const auto & var : vars)
3336 tag_Rn[var->number()],
3337 var->dofIndicesNeighbor(),
3338 var->arrayScalingFactor(),
3345 for (
const auto & vector_tag : vector_tags)
3354 "Non-residual tag in Assembly::addResidualLower");
3359 for (
const auto & var : vars)
3361 tag_Rl[var->number()],
3362 var->dofIndicesLower(),
3363 var->arrayScalingFactor(),
3370 for (
const auto & vector_tag : vector_tags)
3380 "Non-residual tag in Assembly::addResidualScalar");
3386 for (
const auto & var : vars)
3388 residual, tag_Re[var->number()], var->dofIndices(), var->arrayScalingFactor(),
false);
3394 for (
const auto & vector_tag : vector_tags)
3403 for (
const auto & var : vars)
3404 for (
const auto & vector_tag : tags)
3408 _sub_Re[vector_tag._type_id][var->number()],
3410 var->arrayScalingFactor(),
3428 for (
auto & tag : tags)
3434 const std::vector<dof_id_type> & dof_index,
3443 for (MooseIndex(dof_index) i = 0; i < dof_index.size(); ++i)
3454 for (
const auto & var : vars)
3455 for (
const auto & vector_tag : tags)
3459 _sub_Rn[vector_tag._type_id][var->number()],
3460 var->dofIndicesNeighbor(),
3461 var->arrayScalingFactor(),
3469 for (
const auto & var : vars)
3470 for (
const auto & vector_tag : tags)
3474 _sub_Rl[vector_tag._type_id][var->number()],
3475 var->dofIndicesLower(),
3476 var->arrayScalingFactor(),
3483 for (
const auto & vector_tag : tags)
3509 mooseAssert(values.size() == rows.size(),
3510 "Number of cached residuals and number of rows must match!");
3533 mooseAssert(values.size() == rows.size(),
3534 "Number of cached residuals and number of rows must match!");
3536 if (!values.empty())
3548 for (
const auto & var : vars)
3550 tag_Re[var->number()],
3552 var->arrayScalingFactor(),
3563 for (
const auto & var : vars)
3565 tag_Rn[var->number()],
3566 var->dofIndicesNeighbor(),
3567 var->arrayScalingFactor(),
3577 const std::vector<dof_id_type> & idof_indices,
3578 const std::vector<dof_id_type> & jdof_indices)
3580 if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3582 if (jac_block.
n() == 0 || jac_block.
m() == 0)
3587 for (
unsigned int i = 0; i < ivar.
count(); ++i)
3589 unsigned int iv = ivar.
number();
3590 for (
const auto & jt : ConstCouplingRow(iv + i, *
_cm))
3592 unsigned int jv = jvar.
number();
3593 if (jt < jv || jt >= jv + jvar.
count())
3595 unsigned int j = jt - jv;
3599 auto indof = di.size();
3600 auto jndof = dj.size();
3602 unsigned int jj = j;
3607 auto sub = jac_block.
sub_matrix(i * indof, indof, jj * jndof, jndof);
3608 if (scaling_factor[i] != 1.0)
3609 sub *= scaling_factor[i];
3615 _dof_map.constrain_element_matrix(sub, di, dj,
false);
3627 const std::vector<dof_id_type> & idof_indices,
3628 const std::vector<dof_id_type> & jdof_indices,
3631 if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3633 if (jac_block.
n() == 0 || jac_block.
m() == 0)
3640 for (
unsigned int i = 0; i < ivar.
count(); ++i)
3642 unsigned int iv = ivar.
number();
3643 for (
const auto & jt : ConstCouplingRow(iv + i, *
_cm))
3645 unsigned int jv = jvar.
number();
3646 if (jt < jv || jt >= jv + jvar.
count())
3648 unsigned int j = jt - jv;
3652 auto indof = di.size();
3653 auto jndof = dj.size();
3655 unsigned int jj = j;
3660 auto sub = jac_block.
sub_matrix(i * indof, indof, jj * jndof, jndof);
3661 if (scaling_factor[i] != 1.0)
3662 sub *= scaling_factor[i];
3668 _dof_map.constrain_element_matrix(sub, di, dj,
false);
3670 for (MooseIndex(di) i = 0; i < di.size(); i++)
3671 for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3688 const std::vector<dof_id_type> & idof_indices,
3689 const std::vector<dof_id_type> & jdof_indices,
3692 if (idof_indices.size() == 0 || jdof_indices.size() == 0)
3694 if (jac_block.
n() == 0 || jac_block.
m() == 0)
3701 for (
unsigned int i = 0; i < ivar.
count(); ++i)
3703 unsigned int iv = ivar.
number();
3704 for (
const auto & jt : ConstCouplingRow(iv + i, *
_cm))
3706 unsigned int jv = jvar.
number();
3707 if (jt < jv || jt >= jv + jvar.
count())
3709 unsigned int j = jt - jv;
3713 auto indof = di.size();
3714 auto jndof = dj.size();
3716 unsigned int jj = j;
3721 auto sub = jac_block.
sub_matrix(i * indof, indof, jj * jndof, jndof);
3722 if (scaling_factor[i] != 1.0)
3723 sub *= scaling_factor[i];
3725 _dof_map.constrain_element_matrix(sub, di, dj,
false);
3727 for (MooseIndex(di) i = 0; i < di.size(); i++)
3728 for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3729 if (sub(i, j) != 0.0)
3744 const std::vector<dof_id_type> & idof_indices,
3745 const std::vector<dof_id_type> & jdof_indices,
3746 Real scaling_factor,
3751 if ((idof_indices.size() > 0) && (jdof_indices.size() > 0) && jac_block.
n() && jac_block.
m() &&
3754 std::vector<dof_id_type> di(idof_indices);
3755 std::vector<dof_id_type> dj(jdof_indices);
3761 _dof_map.constrain_element_matrix(jac_block, di, dj,
false);
3763 if (scaling_factor != 1.0)
3764 jac_block *= scaling_factor;
3766 for (MooseIndex(di) i = 0; i < di.size(); i++)
3767 for (MooseIndex(dj) j = 0; j < dj.size(); j++)
3781 std::unique_ptr<FEBase> fe(FEBase::build(
elem->dim(), fe_type));
3784 const std::vector<Real> &
JxW = fe->get_JxW();
3785 const std::vector<Point> & q_points = fe->get_xyz();
3789 QGauss qrule(
elem->dim(), fe_type.default_quadrature_order());
3790 fe->attach_quadrature_rule(&qrule);
3795 mooseAssert(qrule.n_points() == q_points.size(),
3796 "The number of points in the quadrature rule doesn't match the number of passed-in " 3797 "points in Assembly::setCoordinateTransformation");
3801 for (
unsigned int qp = 0; qp < qrule.n_points(); ++qp)
3805 vol +=
JxW[qp] * coord;
3816 "Error: Cached data sizes MUST be the same!");
3819 "Error: Cached data sizes MUST be the same for a given tag!");
3883 auto ivar = it.first;
3884 auto jvar = it.second;
3885 auto i = ivar->number();
3886 auto j = jvar->number();
3896 jvar->allDofIndices());
3905 auto ivar = it.first;
3906 auto jvar = it.second;
3907 auto i = ivar->number();
3908 auto j = jvar->number();
3919 jvar->dofIndicesNeighbor());
3925 ivar->dofIndicesNeighbor(),
3926 jvar->dofIndices());
3932 ivar->dofIndicesNeighbor(),
3933 jvar->dofIndicesNeighbor());
3943 auto ivar = it.first;
3944 auto jvar = it.second;
3945 auto i = ivar->number();
3946 auto j = jvar->number();
3955 ivar->dofIndicesLower(),
3956 jvar->dofIndicesLower());
3962 ivar->dofIndicesLower(),
3963 jvar->dofIndicesNeighbor());
3969 ivar->dofIndicesLower(),
3970 jvar->dofIndices());
3976 ivar->dofIndicesNeighbor(),
3977 jvar->dofIndicesLower());
3984 jvar->dofIndicesLower());
3997 jvar->dofIndicesNeighbor());
4003 ivar->dofIndicesNeighbor(),
4004 jvar->dofIndices());
4010 ivar->dofIndicesNeighbor(),
4011 jvar->dofIndicesNeighbor());
4021 auto ivar = it.first;
4022 auto jvar = it.second;
4023 auto i = ivar->number();
4024 auto j = jvar->number();
4033 ivar->dofIndicesLower(),
4034 jvar->dofIndicesLower());
4040 ivar->dofIndicesLower(),
4041 jvar->dofIndices());
4048 jvar->dofIndicesLower());
4088 auto ivar = it.first;
4089 auto jvar = it.second;
4090 auto i = ivar->number();
4091 auto j = jvar->number();
4100 jvar->allDofIndices(),
4110 auto ivar = it.first;
4111 auto jvar = it.second;
4112 auto i = ivar->number();
4113 auto j = jvar->number();
4124 jvar->dofIndicesNeighbor(),
4129 ivar->dofIndicesNeighbor(),
4136 ivar->dofIndicesNeighbor(),
4137 jvar->dofIndicesNeighbor(),
4148 auto ivar = it.first;
4149 auto jvar = it.second;
4150 auto i = ivar->number();
4151 auto j = jvar->number();
4159 ivar->dofIndicesLower(),
4160 jvar->dofIndicesLower(),
4166 ivar->dofIndicesLower(),
4173 ivar->dofIndicesLower(),
4174 jvar->dofIndicesNeighbor(),
4181 jvar->dofIndicesLower(),
4196 jvar->dofIndicesNeighbor(),
4202 ivar->dofIndicesNeighbor(),
4203 jvar->dofIndicesLower(),
4209 ivar->dofIndicesNeighbor(),
4216 ivar->dofIndicesNeighbor(),
4217 jvar->dofIndicesNeighbor(),
4227 const DofMap & dof_map,
4228 std::vector<dof_id_type> & dof_indices,
4230 const std::set<TagID> & tags)
4232 for (
auto tag : tags)
4240 const DofMap & dof_map,
4241 std::vector<dof_id_type> & dof_indices,
4245 if (dof_indices.size() == 0)
4247 if (!(*
_cm)(ivar, jvar))
4254 const unsigned int ivn = iv.number();
4255 const unsigned int jvn = jv.number();
4264 const unsigned int i = ivar - ivn;
4265 const unsigned int j = jvar - jvn;
4268 auto di = dof_indices;
4269 auto dj = dof_indices;
4271 auto indof = di.size();
4272 auto jndof = dj.size();
4274 unsigned int jj = j;
4278 auto sub = ke.
sub_matrix(i * indof, indof, jj * jndof, jndof);
4283 dof_map.constrain_element_matrix(sub, di, dj,
false);
4285 if (scaling_factor[i] != 1.0)
4286 sub *= scaling_factor[i];
4293 const unsigned int ivar,
4294 const unsigned int jvar,
4295 const DofMap & dof_map,
4296 const std::vector<dof_id_type> & idof_indices,
4297 const std::vector<dof_id_type> & jdof_indices,
4301 if (idof_indices.size() == 0 || jdof_indices.size() == 0)
4303 if (jacobian.
n() == 0 || jacobian.
m() == 0)
4305 if (!(*
_cm)(ivar, jvar))
4312 const unsigned int ivn = iv.number();
4313 const unsigned int jvn = jv.number();
4322 const unsigned int i = ivar - ivn;
4323 const unsigned int j = jvar - jvn;
4326 auto di = idof_indices;
4327 auto dj = jdof_indices;
4329 auto indof = di.size();
4330 auto jndof = dj.size();
4332 unsigned int jj = j;
4336 auto sub = keg.
sub_matrix(i * indof, indof, jj * jndof, jndof);
4341 dof_map.constrain_element_matrix(sub, di, dj,
false);
4343 if (scaling_factor[i] != 1.0)
4344 sub *= scaling_factor[i];
4351 const unsigned int ivar,
4352 const unsigned int jvar,
4353 const DofMap & dof_map,
4354 const std::vector<dof_id_type> & idof_indices,
4355 const std::vector<dof_id_type> & jdof_indices,
4357 const std::set<TagID> & tags)
4359 for (
auto tag : tags)
4361 jacobian, ivar, jvar, dof_map, idof_indices, jdof_indices,
GlobalDataKey{}, tag);
4366 const unsigned int ivar,
4367 const unsigned int jvar,
4368 const DofMap & dof_map,
4369 std::vector<dof_id_type> & dof_indices,
4370 std::vector<dof_id_type> & neighbor_dof_indices,
4374 if (dof_indices.size() == 0 && neighbor_dof_indices.size() == 0)
4376 if (!(*
_cm)(ivar, jvar))
4383 const unsigned int ivn = iv.number();
4384 const unsigned int jvn = jv.number();
4395 const unsigned int i = ivar - ivn;
4396 const unsigned int j = jvar - jvn;
4398 auto dc = dof_indices;
4399 auto dn = neighbor_dof_indices;
4400 auto cndof = dc.size();
4401 auto nndof = dn.size();
4403 unsigned int jj = j;
4407 auto suben = ken.
sub_matrix(i * cndof, cndof, jj * nndof, nndof);
4408 auto subne = kne.
sub_matrix(i * nndof, nndof, jj * cndof, cndof);
4409 auto subnn = knn.
sub_matrix(i * nndof, nndof, jj * nndof, nndof);
4416 dof_map.constrain_element_matrix(suben, dc, dn,
false);
4417 dof_map.constrain_element_matrix(subne, dn, dc,
false);
4418 dof_map.constrain_element_matrix(subnn, dn, dn,
false);
4421 if (scaling_factor[i] != 1.0)
4423 suben *= scaling_factor[i];
4424 subne *= scaling_factor[i];
4425 subnn *= scaling_factor[i];
4435 const unsigned int ivar,
4436 const unsigned int jvar,
4437 const DofMap & dof_map,
4438 std::vector<dof_id_type> & dof_indices,
4439 std::vector<dof_id_type> & neighbor_dof_indices,
4441 const std::set<TagID> & tags)
4443 for (
const auto tag : tags)
4445 jacobian, ivar, jvar, dof_map, dof_indices, neighbor_dof_indices,
GlobalDataKey{}, tag);
4460 for (
const auto & var_j : vars)
4478 const std::set<TagID> & tags)
4480 for (
auto tag : tags)
4529 mooseAssert(
_xfem !=
nullptr,
"This function should not be called if xfem is inactive");
4538 "Size of weight multipliers in xfem doesn't match number of quadrature points");
4539 for (
unsigned i = 0; i < xfem_weight_multipliers.
size(); i++)
4542 xfem_weight_multipliers.
release();
4549 mooseAssert(
_xfem !=
nullptr,
"This function should not be called if xfem is inactive");
4555 if (
_xfem->getXFEMFaceWeights(
4559 "Size of weight multipliers in xfem doesn't match number of quadrature points");
4560 for (
unsigned i = 0; i < xfem_face_weight_multipliers.
size(); i++)
4563 xfem_face_weight_multipliers.
release();
4579 for (MooseIndex(weights.size()) i = 0; i < weights.size(); ++i)
4585 Assembly::fePhi<VectorValue<Real>>(FEType type)
const 4587 buildVectorFE(type);
4588 return _vector_fe_shape_data[type]->_phi;
4593 Assembly::feGradPhi<VectorValue<Real>>(FEType type)
const 4595 buildVectorFE(type);
4596 return _vector_fe_shape_data[type]->_grad_phi;
4601 Assembly::feSecondPhi<VectorValue<Real>>(FEType type)
const 4603 _need_second_derivative[type] =
true;
4604 buildVectorFE(type);
4605 return _vector_fe_shape_data[type]->_second_phi;
4610 Assembly::fePhiLower<VectorValue<Real>>(FEType type)
const 4612 buildVectorLowerDFE(type);
4613 return _vector_fe_shape_data_lower[type]->_phi;
4618 Assembly::feDualPhiLower<VectorValue<Real>>(FEType type)
const 4620 buildVectorDualLowerDFE(type);
4621 return _vector_fe_shape_data_dual_lower[type]->_phi;
4626 Assembly::feGradPhiLower<VectorValue<Real>>(FEType type)
const 4628 buildVectorLowerDFE(type);
4629 return _vector_fe_shape_data_lower[type]->_grad_phi;
4634 Assembly::feGradDualPhiLower<VectorValue<Real>>(FEType type)
const 4636 buildVectorDualLowerDFE(type);
4637 return _vector_fe_shape_data_dual_lower[type]->_grad_phi;
4642 Assembly::fePhiFace<VectorValue<Real>>(FEType type)
const 4644 buildVectorFaceFE(type);
4645 return _vector_fe_shape_data_face[type]->_phi;
4650 Assembly::feGradPhiFace<VectorValue<Real>>(FEType type)
const 4652 buildVectorFaceFE(type);
4653 return _vector_fe_shape_data_face[type]->_grad_phi;
4658 Assembly::feSecondPhiFace<VectorValue<Real>>(FEType type)
const 4660 _need_second_derivative[type] =
true;
4661 buildVectorFaceFE(type);
4662 return _vector_fe_shape_data_face[type]->_second_phi;
4667 Assembly::fePhiNeighbor<VectorValue<Real>>(FEType type)
const 4669 buildVectorNeighborFE(type);
4670 return _vector_fe_shape_data_neighbor[type]->_phi;
4675 Assembly::feGradPhiNeighbor<VectorValue<Real>>(FEType type)
const 4677 buildVectorNeighborFE(type);
4678 return _vector_fe_shape_data_neighbor[type]->_grad_phi;
4683 Assembly::feSecondPhiNeighbor<VectorValue<Real>>(FEType type)
const 4685 _need_second_derivative_neighbor[type] =
true;
4686 buildVectorNeighborFE(type);
4687 return _vector_fe_shape_data_neighbor[type]->_second_phi;
4692 Assembly::fePhiFaceNeighbor<VectorValue<Real>>(FEType type)
const 4694 buildVectorFaceNeighborFE(type);
4695 return _vector_fe_shape_data_face_neighbor[type]->_phi;
4700 Assembly::feGradPhiFaceNeighbor<VectorValue<Real>>(FEType type)
const 4702 buildVectorFaceNeighborFE(type);
4703 return _vector_fe_shape_data_face_neighbor[type]->_grad_phi;
4708 Assembly::feSecondPhiFaceNeighbor<VectorValue<Real>>(FEType type)
const 4710 _need_second_derivative_neighbor[type] =
true;
4711 buildVectorFaceNeighborFE(type);
4712 return _vector_fe_shape_data_face_neighbor[type]->_second_phi;
4717 Assembly::feCurlPhi<VectorValue<Real>>(FEType type)
const 4719 _need_curl[type] =
true;
4720 buildVectorFE(type);
4721 return _vector_fe_shape_data[type]->_curl_phi;
4726 Assembly::feCurlPhiFace<VectorValue<Real>>(FEType type)
const 4728 _need_curl[type] =
true;
4729 buildVectorFaceFE(type);
4730 return _vector_fe_shape_data_face[type]->_curl_phi;
4735 Assembly::feCurlPhiNeighbor<VectorValue<Real>>(FEType type)
const 4737 _need_curl[type] =
true;
4738 buildVectorNeighborFE(type);
4739 return _vector_fe_shape_data_neighbor[type]->_curl_phi;
4744 Assembly::feCurlPhiFaceNeighbor<VectorValue<Real>>(FEType type)
const 4746 _need_curl[type] =
true;
4747 buildVectorFaceNeighborFE(type);
4748 return _vector_fe_shape_data_face_neighbor[type]->_curl_phi;
4753 Assembly::feDivPhi<VectorValue<Real>>(FEType type)
const 4755 _need_div[type] =
true;
4756 buildVectorFE(type);
4757 return _vector_fe_shape_data[type]->_div_phi;
4762 Assembly::feDivPhiFace<VectorValue<Real>>(FEType type)
const 4764 _need_div[type] =
true;
4765 buildVectorFaceFE(type);
4766 return _vector_fe_shape_data_face[type]->_div_phi;
4771 Assembly::feDivPhiNeighbor<VectorValue<Real>>(FEType type)
const 4773 _need_div[type] =
true;
4774 buildVectorNeighborFE(type);
4775 return _vector_fe_shape_data_neighbor[type]->_div_phi;
4780 Assembly::feDivPhiFaceNeighbor<VectorValue<Real>>(FEType type)
const 4782 _need_div[type] =
true;
4783 buildVectorFaceNeighborFE(type);
4784 return _vector_fe_shape_data_face_neighbor[type]->_div_phi;
4792 const FEType helper_type(helper_order,
LAGRANGE);
4795 feSecondPhi<Real>(helper_type);
4796 feSecondPhiFace<Real>(helper_type);
4840 const std::unordered_set<FEFamily> disable_families(disable_p_refinement_for_families.begin(),
4841 disable_p_refinement_for_families.end());
4844 const FEType helper_type(helper_order,
LAGRANGE);
4846 [&disable_families](
const unsigned int num_dimensionalities,
auto & fe_container)
4848 if (!disable_families.empty())
4851 auto fe_container_it = fe_container.find(
dim);
4852 if (fe_container_it != fe_container.end())
4853 for (
auto & [fe_type, fe_ptr] : fe_container_it->second)
4854 if (disable_families.count(fe_type.family))
4855 fe_ptr->add_p_level_in_reinit(
false);
4858 auto process_fe_and_helpers = [process_fe, &helper_type](
auto & unique_helper_container,
4859 auto & helper_container,
4860 const unsigned int num_dimensionalities,
4861 const bool user_added_helper_type,
4862 auto & fe_container)
4864 unique_helper_container.resize(num_dimensionalities);
4867 auto & unique_helper = unique_helper_container[
dim];
4868 unique_helper = FEGenericBase<Real>::build(
dim, helper_type);
4870 unique_helper->add_p_level_in_reinit(
false);
4871 helper_container[
dim] = unique_helper.get();
4876 if (!user_added_helper_type)
4878 auto & fe_container_dim = libmesh_map_find(fe_container,
dim);
4879 auto fe_it = fe_container_dim.find(helper_type);
4880 mooseAssert(fe_it != fe_container_dim.end(),
"We should have the helper type");
4881 delete fe_it->second;
4882 fe_container_dim.erase(fe_it);
4886 process_fe(num_dimensionalities, fe_container);
4929 const Point & point,
4939 const Point & point,
4950 Assembly::genericQPoints<false>()
const 4957 Assembly::genericQPoints<true>()
const const Elem *const & elem() const
Return the current element.
virtual MooseMesh & mesh()=0
MooseArray< DualReal > _ad_JxW
void copyShapes(MooseVariableField< T > &v)
void setWeights(const std::vector< Real > &weights)
Set the quadrature weights.
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face
types of finite elements
bool _need_neighbor_elem_volume
true is apps need to compute neighbor element volume
void cacheResidualNeighbor(GlobalDataKey, const std::vector< VectorTag > &tags)
Takes the values that are currently in _sub_Rn of all field variables and appends them to the cached ...
virtual void insert(const Number *v, const std::vector< numeric_index_type > &dof_indices)
ArbitraryQuadrature * _current_qrule_arbitrary
The current arbitrary quadrature rule used within the element interior.
std::map< FEType, FEBase * > _current_fe
The "volume" fe object that matches the current elem.
MooseArray< Real > _curvatures
const std::vector< MooseVariableFieldBase * > & getVariables(THREAD_ID tid)
VectorVariablePhiValue _phi
std::vector< std::vector< dof_id_type > > _cached_residual_rows
Where the cached values should go (the first vector is for TIME vs NONTIME)
std::vector< std::vector< dof_id_type > > _cached_jacobian_rows
Row where the corresponding cached value should go.
void computeGradPhiAD(const Elem *elem, unsigned int n_qp, ADTemplateVariablePhiGradient< OutputType > &grad_phi, FEGenericBase< OutputType > *fe)
compute gradient of phi possibly with derivative information with respect to nonlinear displacement v...
unsigned int _max_cached_residuals
std::map< FEType, ADTemplateVariablePhiGradient< Real > > _ad_grad_phi_data_face
std::vector< dof_id_type > componentDofIndices(const std::vector< dof_id_type > &dof_indices, unsigned int component) const
Obtain DoF indices of a component with the indices of the 0th component.
MooseArray< VectorValue< DualReal > > _ad_q_points_face
std::vector< DualReal > _ad_detadz_map
QBase * _current_qrule_lower
quadrature rule used on lower dimensional elements.
bool _user_added_fe_lower_of_helper_type
MooseArray< Point > _current_physical_points
This will be filled up with the physical points passed into reinitAtPhysical() if it is called...
std::vector< DualReal > _ad_dxidy_map
const VariablePhiGradient & gradPhiFaceNeighbor(const MooseVariableField< Real > &) const
std::vector< std::unique_ptr< FEBase > > _unique_fe_lower_helper
MooseArray< VectorValue< DualReal > > _ad_normals
virtual const FieldVariablePhiValue & phi() const =0
Return the variable's elemental shape functions.
void setResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Set a local residual block to a global residual vector with proper scaling.
void buildNeighborFE(FEType type) const
Build FEs for a neighbor with a type.
const VariablePhiSecond & secondPhi() const
void setMortarQRule(Order order)
Specifies a custom qrule for integration on mortar segment mesh.
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe
Each dimension's actual vector fe objects indexed on type.
typename OutputTools< typename Moose::ADType< T >::type >::VariablePhiGradient ADTemplateVariablePhiGradient
void buildVectorLowerDFE(FEType type) const
Build Vector FEs for a lower dimensional element with a type.
std::vector< VectorValue< DualReal > > _ad_d2xyzdxi2_map
MooseArray< Real > _coord_neighbor
The current coordinate transformation coefficients.
virtual void haveADObjects(bool have_ad_objects)
Method for setting whether we have any ad objects.
void buildFE(FEType type) const
Build FEs with a type.
QBase * _qrule_msm
A qrule object for working on mortar segement elements.
const std::vector< MooseVariableScalar * > & getScalarVariables(THREAD_ID tid)
QBase * _current_qrule
The current current quadrature rule being used (could be either volumetric or arbitrary - for dirac k...
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kle
dlower/dsecondary (or dlower/delement)
virtual void zero() override final
void cacheResidualLower(GlobalDataKey, const std::vector< VectorTag > &tags)
Takes the values that are currently in _sub_Rl and appends them to the cached values.
const unsigned int invalid_uint
bool hasVector(const std::string &tag_name) const
Check if the named vector exists in the system.
MooseArray< DualReal > _ad_JxW_face
Assembly(SystemBase &sys, THREAD_ID tid)
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data
Shape function values, gradients, second derivatives for each FE type.
void buildLowerDDualFE(FEType type) const
void addJacobianBlockNonlocal(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, GlobalDataKey, TagID tag)
Adds non-local element matrix for ivar rows and jvar columns to the global Jacobian matrix...
MooseArray< Real > _current_JxW_neighbor
The current transformed jacobian weights on a neighbor's face.
std::shared_ptr< XFEMInterface > _xfem
The XFEM controller.
void reinitNeighborLowerDElem(const Elem *elem)
reinitialize a neighboring lower dimensional element
void prepareJacobianBlock()
Sizes and zeroes the Jacobian blocks used for the current element.
DenseMatrix< Number > & jacobianBlockNonlocal(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Get local Jacobian block from non-local contribution for a pair of variables and a tag...
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
virtual void zero() override final
std::vector< DualReal > _ad_dzetady_map
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kel
dsecondary/dlower (or delement/dlower)
const std::vector< Real > & arrayScalingFactor() const
MooseArray< Real > _coord
The current coordinate transformation coefficients.
unsigned int number() const
Get variable number coming from libMesh.
MooseArray< VectorValue< DualReal > > _ad_q_points
std::unique_ptr< QBase > vol
volume/elem (meshdim) quadrature rule
DenseMatrix< Number > & jacobianBlockNeighbor(Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Get local Jacobian block of a DG Jacobian type for a pair of variables and a tag. ...
void reinitFE(const Elem *elem)
Just an internal helper function to reinit the volume FE objects.
TagID _id
The id associated with the vector tag.
virtual unsigned int currentNlSysNum() const =0
void jacobianBlockUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not Jacobian coupling between ivar and jvar is used to the value used.
void coordTransformFactor(const P &point, C &factor, const Moose::CoordinateSystemType coord_type, const unsigned int rz_radial_coord=libMesh::invalid_uint)
compute a coordinate transformation factor
void init(const CouplingMatrix *cm)
Initialize the Assembly object and set the CouplingMatrix for use throughout.
std::map< unsigned int, FEBase * > _holder_fe_neighbor_helper
Each dimension's helper objects.
std::unique_ptr< FEBase > _fe_msm
A FE object for working on mortar segement elements.
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_face_neighbor
void createQRules(QuadratureType type, Order order, Order volume_order, Order face_order, SubdomainID block, bool allow_negative_qweights=true)
Creates block-specific volume, face and arbitrary qrules based on the orders and the flag of whether ...
std::vector< DualReal > _ad_dzetadz_map
const Elem & elem() const
void reinitAtPhysical(const Elem *elem, const std::vector< Point > &physical_points)
Reinitialize the assembly data at specific physical point in the given element.
void addCachedResidualDirectly(NumericVector< Number > &residual, GlobalDataKey, const VectorTag &vector_tag)
Adds the values that have been cached by calling cacheResidual(), cacheResidualNeighbor(), and/or cacheResidualLower() to a user-defined residual (that is, not necessarily the vector that vector_tag points to)
bool _current_elem_volume_computed
Boolean to indicate whether current element volumes has been computed.
void jacobianBlockLowerUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not lower Jacobian coupling between ivar and jvar is used to the value used...
const VectorVariablePhiDivergence & divPhi(const MooseVariableField< RealVectorValue > &) const
std::map< FEType, bool > _need_div
std::vector< std::unique_ptr< FEBase > > _unique_fe_neighbor_helper
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kll
dlower/dlower
bool _user_added_fe_face_neighbor_of_helper_type
void setResidualNeighbor(NumericVector< Number > &residual, GlobalDataKey, const VectorTag &vector_tag)
Sets local neighbor residuals of all field variables to the global residual vector for a tag...
bool _have_p_refinement
Whether we have ever conducted p-refinement.
Real _current_neighbor_volume
Volume of the current neighbor.
void jacobianBlockNonlocalUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not nonlocal Jacobian coupling between ivar and jvar is used to the value used...
virtual void add_vector(const Number *v, const std::vector< numeric_index_type > &dof_indices)
ArbitraryQuadrature * qruleArbitraryFace(const Elem *elem, unsigned int side)
std::map< FEType, FEBase * > _current_fe_face
The "face" fe object that matches the current elem.
const VectorVariablePhiCurl & curlPhi(const MooseVariableField< RealVectorValue > &) const
std::vector< Point > _current_neighbor_ref_points
The current reference points on the neighbor element.
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kln
dlower/dprimary (or dlower/dneighbor)
const Elem * _current_neighbor_elem
The current neighbor "element".
void addJacobianNonlocal(GlobalDataKey)
Adds non-local Jacobian to the global Jacobian matrices.
std::vector< std::vector< DenseVector< Number > > > _sub_Rn
unsigned int count() const
Get the number of components Note: For standard and vector variables, the number is one...
virtual const FieldVariablePhiGradient & gradPhiNeighbor() const =0
Return the gradients of the variable's shape functions on a neighboring element.
const VariablePhiValue & phi() const
TagTypeID _type_id
The index for this tag into a vector that contains tags of only its type ordered by ID...
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_neighbor
std::map< unsigned int, FEBase * > _holder_fe_face_neighbor_helper
void setPoints(const std::vector< Point > &points)
Set the quadrature points.
const VariablePhiGradient & gradPhi() const
MooseArray< Real > _current_JxW_face
The current transformed jacobian weights on a face.
void modifyWeightsDueToXFEM(const Elem *elem)
Update the integration weights for XFEM partial elements.
void coordTransformFactorRZGeneral(const P &point, const std::pair< Point, RealVectorValue > &axis, C &factor)
Computes a coordinate transformation factor for a general axisymmetric axis.
const Elem * _current_elem
The current "element" we are currently on.
const VariablePhiSecond & secondPhiNeighbor(const MooseVariableField< Real > &) const
bool computingScalingJacobian() const
Whether we are computing an initial Jacobian for automatic variable scaling.
THREAD_ID _tid
Thread number (id)
void reinitFEFaceNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
virtual const FieldVariablePhiSecond & secondPhi() const =0
Return the rank-2 tensor of second derivatives of the variable's elemental shape functions.
void reinitNeighborAtPhysical(const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > &physical_points)
Reinitializes the neighbor at the physical coordinates on neighbor side given.
std::vector< std::unique_ptr< FEBase > > _unique_fe_face_helper
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_nonlocal_used
std::vector< std::pair< MooseVariableScalar *, MooseVariableFieldBase * > > _cm_sf_entry
Entries in the coupling matrix for scalar variables vs field variables.
void resizeADMappingObjects(unsigned int n_qp, unsigned int dim)
resize any objects that contribute to automatic differentiation-related mapping calculations ...
OutputTools< Real >::VariablePhiValue VariablePhiValue
virtual bool hasMatrix(TagID tag) const
Check if the tagged matrix exists in the system.
ElemSideBuilder _current_neighbor_side_elem_builder
In place side element builder for _current_neighbor_side_elem.
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_face_neighbor
std::map< FEType, bool > _need_second_derivative
const VariablePhiValue & phiFaceNeighbor(const MooseVariableField< Real > &) const
unsigned int elemSideID() const
std::vector< bool > _component_block_diagonal
An flag array Indiced by variable index to show if there is no component-wise coupling for the variab...
Real _current_elem_volume
Volume of the current element.
std::map< FEType, FEBase * > _current_fe_face_neighbor
The "neighbor face" fe object that matches the current elem.
const MooseArray< ADReal > & adCurvatures() const
This class provides an interface for common operations on field variables of both FE and FV types wit...
std::map< FEType, FEVectorBase * > _current_vector_fe_face
The "face" vector fe object that matches the current elem.
void shallowCopy(const MooseArray &rhs)
Doesn't actually make a copy of the data.
std::vector< std::unique_ptr< FEBase > > _unique_fe_helper
Containers for holding unique FE helper types if we are doing p-refinement.
std::vector< DualReal > _ad_dzetadx_map
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
void processLocalResidual(DenseVector< Number > &res_block, std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Appling scaling, constraints to the local residual block and populate the full DoF indices for array ...
const FEType _helper_type
The finite element type of the FE helper classes.
void addResidualLower(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add local neighbor residuals of all field variables for a set of tags onto the global residual vector...
DenseMatrix< Number > & jacobianBlock(unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Get local Jacobian block for a pair of variables and a tag.
void addJacobianScalar(GlobalDataKey)
Add Jacobians for pairs of scalar variables into the global Jacobian matrices.
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Knn
jacobian contributions from the neighbor <Tag, ivar, jvar>
Base class for a system (of equations)
const VariablePhiValue & phiNeighbor(const MooseVariableField< Real > &) const
void setVolumeQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for volume integration.
FEGenericBase< RealGradient > FEVectorBase
unsigned int numExtraElemIntegers() const
Number of extra element integers Assembly tracked.
std::vector< std::vector< Real > > _cached_jacobian_values
Values cached by calling cacheJacobian()
std::vector< std::pair< MooseVariableFieldBase *, MooseVariableFieldBase * > > _cm_ff_entry
Entries in the coupling matrix for field variables.
void addJacobianNeighborLowerD(GlobalDataKey)
Add all portions of the Jacobian except PrimaryPrimary, e.g.
std::map< FEType, FEVectorBase * > _current_vector_fe
The "volume" vector fe object that matches the current elem.
VariablePhiGradient _grad_phi
template void coordTransformFactor< ADPoint, ADReal >(const SubProblem &s, SubdomainID sub_id, const ADPoint &point, ADReal &factor, SubdomainID neighbor_sub_id)
virtual const FieldVariablePhiValue & phiNeighbor() const =0
Return the variable's shape functions on a neighboring element.
ElemSideBuilder _compute_face_map_side_elem_builder
In place side element builder for computeFaceMap()
void modifyFaceWeightsDueToXFEM(const Elem *elem, unsigned int side=0)
Update the face integration weights for XFEM partial elements.
const CouplingMatrix * _cm
Coupling matrices.
std::map< FEType, bool > _need_second_derivative_neighbor
void reinitElemAndNeighbor(const Elem *elem, unsigned int side, const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > *neighbor_reference_points=nullptr)
Reinitialize an element and its neighbor along a particular side.
OutputTools< Real >::VariablePhiSecond VariablePhiSecond
unsigned int neighborSideID() const
DenseVector< Number > _tmp_Re
auxiliary vector for scaling residuals (optimization to avoid expensive construction/destruction) ...
unsigned int _mesh_dimension
std::vector< Eigen::Map< RealDIMValue > > _mapped_normals
Mapped normals.
VectorVariablePhiDivergence _vector_div_phi_face
void cacheJacobian(GlobalDataKey)
Takes the values that are currently in _sub_Kee and appends them to the cached values.
void reinit(const Elem *elem)
Reinitialize objects (JxW, q_points, ...) for an elements.
const NumericVector< Real > * _scaling_vector
The map from global index to variable scaling factor.
auto max(const L &left, const R &right)
QBase * _current_qrule_volume
The current volumetric quadrature for the element.
virtual void set(const numeric_index_type i, const numeric_index_type j, const Number value)=0
void addCachedResiduals(GlobalDataKey, const std::vector< VectorTag > &tags)
Pushes all cached residuals to the global residual vectors associated with each tag.
const DofMap & _dof_map
DOF map.
Data structure for tracking/grouping a set of quadrature rules for a particular dimensionality of mes...
std::vector< DualReal > _ad_dxidx_map
unsigned int _max_cached_jacobians
virtual void add(const numeric_index_type i, const numeric_index_type j, const Number value)=0
DenseMatrix< Number > & jacobianBlockMortar(Moose::ConstraintJacobianType type, unsigned int ivar, unsigned int jvar, LocalDataKey, TagID tag)
Returns the jacobian block for the given mortar Jacobian type.
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Ken
jacobian contributions from the element and neighbor <Tag, ivar, jvar>
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_dual_lower
virtual unsigned int nVariables() const
Get the number of variables in this system.
const CouplingMatrix & _nonlocal_cm
virtual const FieldVariablePhiValue & phiFaceNeighbor() const =0
Return the variable's shape functions on a neighboring element face.
std::vector< Point > _temp_reference_points
Temporary work data for reinitAtPhysical()
std::vector< std::pair< MooseVariableScalar *, MooseVariableScalar * > > _cm_ss_entry
Entries in the coupling matrix for scalar variables.
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kne
jacobian contributions from the neighbor and element <Tag, ivar, jvar>
const Elem * _current_neighbor_lower_d_elem
The current neighboring lower dimensional element.
unsigned int _current_neighbor_side
The current side of the selected neighboring element (valid only when working with sides) ...
void cacheResidualNodes(const DenseVector< Number > &res, const std::vector< dof_id_type > &dof_index, LocalDataKey, TagID tag)
Lets an external class cache residual at a set of nodes.
void prepareLowerD()
Prepare the Jacobians and residuals for a lower dimensional element.
bool _user_added_fe_neighbor_of_helper_type
void buildFaceNeighborFE(FEType type) const
Build FEs for a neighbor face with a type.
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_face
void reinitLowerDElem(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)
Reinitialize FE data for a lower dimenesional element with a given set of reference points...
void computeCurrentFaceVolume()
std::vector< std::vector< DenseVector< Number > > > _sub_Rl
residual contributions for each variable from the lower dimensional element
This data structure is used to store geometric and variable related metadata about each cell face in ...
QRules & qrules(unsigned int dim)
virtual const FieldVariablePhiSecond & secondPhiFaceNeighbor() const =0
Return the rank-2 tensor of second derivatives of the variable's shape functions on a neighboring ele...
virtual const std::vector< dof_id_type > & dofIndicesNeighbor() const =0
Get neighbor DOF indices for currently selected element.
static const subdomain_id_type invalid_subdomain_id
virtual void add_matrix(const DenseMatrix< Number > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
void prepareOffDiagScalar()
VectorVariablePhiSecond _second_phi
unsigned int size() const
The number of elements that can currently be stored in the array.
Implements a fake quadrature rule where you can specify the locations (in the reference domain) of th...
void addJacobianCoupledVarPair(const MooseVariableBase &ivar, const MooseVariableBase &jvar)
Adds element matrices for ivar rows and jvar columns to the global Jacobian matrices.
void addJacobianLowerD(GlobalDataKey)
Add portions of the Jacobian of LowerLower, LowerSecondary, and SecondaryLower for boundary condition...
virtual const FieldVariablePhiSecond & secondPhiFace() const =0
Return the rank-2 tensor of second derivatives of the variable's shape functions on an element face...
QBase * _current_qrule_neighbor
quadrature rule used on neighbors
std::vector< DualReal > _ad_dxidz_map
const Elem * neighborPtr() const
VectorVariablePhiGradient _grad_phi
ArbitraryQuadrature * _current_qrule_arbitrary_face
The current arbitrary quadrature rule used on the element face.
void modifyArbitraryWeights(const std::vector< Real > &weights)
Modify the weights when using the arbitrary quadrature rule.
std::vector< std::pair< MooseVariableFieldBase *, MooseVariableScalar * > > _cm_fs_entry
Entries in the coupling matrix for field variables vs scalar variables.
const std::vector< Real > * _JxW_msm
A JxW for working on mortar segement elements.
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
void coordTransformFactor(const SubProblem &s, const SubdomainID sub_id, const P &point, C &factor, const SubdomainID neighbor_sub_id)
Computes a conversion multiplier for use when computing integraals for the current coordinate system ...
std::unique_ptr< QBase > face
area/face (meshdim-1) quadrature rule
dof_id_type numeric_index_type
std::map< FEType, ADTemplateVariablePhiGradient< Real > > _ad_grad_phi_data
template void coordTransformFactor< Point, Real >(const SubProblem &s, SubdomainID sub_id, const Point &point, Real &factor, SubdomainID neighbor_sub_id)
std::unordered_map< SubdomainID, std::vector< QRules > > _qrules
Holds quadrature rules for each dimension.
void prepareResidual()
Sizes and zeroes the residual for the current element.
VectorVariablePhiCurl _curl_phi
std::map< FEType, FEBase * > _current_fe_neighbor
The "neighbor" fe object that matches the current elem.
Real elementVolume(const Elem *elem) const
On-demand computation of volume element accounting for RZ/RSpherical.
void cacheJacobianNeighbor(GlobalDataKey)
Takes the values that are currently in the neighbor Dense Matrices and appends them to the cached val...
SubdomainID _current_neighbor_subdomain_id
The current neighbor subdomain ID.
std::vector< std::vector< DenseVector< Number > > > _sub_Re
void reinitNeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
Reinitializes the neighbor side using reference coordinates.
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_used
Flag that indicates if the jacobian block was used.
virtual numeric_index_type m() const =0
void setFaceQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for face integration.
std::vector< DualReal > _ad_jac
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
void buildFaceFE(FEType type) const
Build FEs for a face with a type.
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_neighbor
types of finite elements
void cacheResidual(GlobalDataKey, const std::vector< VectorTag > &tags)
Takes the values that are currently in _sub_Re of all field variables and appends them to the cached ...
MooseArray< std::vector< Point > > _current_tangents
The current tangent vectors at the quadrature points.
virtual const std::vector< dof_id_type > & dofIndices() const
Get local DoF indices.
bool _calculate_curvatures
const std::vector< VectorTag > & _residual_vector_tags
The residual vector tags that Assembly could possibly contribute to.
std::map< FEType, FEVectorBase * > _current_vector_fe_face_neighbor
The "neighbor face" vector fe object that matches the current elem.
std::vector< dof_id_type > _neighbor_extra_elem_ids
Extra element IDs of neighbor.
MooseArray< Real > _current_JxW
The current list of transformed jacobian weights.
virtual void zero_rows(std::vector< numeric_index_type > &rows, Number diag_value=0.0)
void havePRefinement(const std::vector< FEFamily > &disable_p_refinement_for_families)
Indicate that we have p-refinement.
virtual bool checkNonlocalCouplingRequirement()
std::vector< std::pair< unsigned int, unsigned short > > _disp_numbers_and_directions
Container of displacement numbers and directions.
std::vector< dof_id_type > _temp_dof_indices
Temporary work vector to keep from reallocating it.
OutputTools< Real >::VariablePhiCurl VariablePhiCurl
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Kee
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_lower
OStreamProxy err(std::cerr)
const VariablePhiSecond & secondPhiFace(const MooseVariableField< Real > &) const
std::map< FEType, ADTemplateVariablePhiGradient< RealVectorValue > > _ad_vector_grad_phi_data
bool usesPhiNeighbor() const
Whether or not this variable is actually using the shape function value.
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_face
types of vector finite elements
subdomain_id_type SubdomainID
virtual bool computingSecond() const =0
Whether or not this variable is computing any second derivatives.
unsigned int number() const
Gets the number of this system.
std::vector< T > stdVector() const
Extremely inefficient way to produce a std::vector from a MooseArray!
QBase * qruleFace(const Elem *elem, unsigned int side)
This is an abstraction over the internal qrules function.
bool _calculate_ad_coord
Whether to calculate coord with AD.
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_lower
const Elem * _current_lower_d_elem
The current lower dimensional element.
void addResidualBlock(NumericVector< Number > &residual, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Add a local residual block to a global residual vector with proper scaling.
void addResidual(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add local residuals of all field variables for a set of tags onto the global residual vectors associa...
void addJacobianNeighborTags(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, std::vector< dof_id_type > &dof_indices, std::vector< dof_id_type > &neighbor_dof_indices, GlobalDataKey, const std::set< TagID > &tags)
Adds three neighboring element matrices for ivar rows and jvar columns to the global Jacobian matrix...
std::vector< VectorValue< DualReal > > _ad_dxyzdeta_map
const VariablePhiValue & phiFace() const
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_lower
FE objects for lower dimensional elements.
bool hasSecondOrderElements()
check if the mesh has SECOND order elements
FEGenericBase< Real > FEBase
void addJacobianBlockNonlocalTags(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, GlobalDataKey, const std::set< TagID > &tags)
Adds non-local element matrix for ivar rows and jvar columns to the global Jacobian matrix...
OutputTools< Real >::VariablePhiDivergence VariablePhiDivergence
Real _current_lower_d_elem_volume
The current lower dimensional element volume.
std::map< unsigned int, std::map< FEType, FEBase * > > _fe
Each dimension's actual fe objects indexed on type.
const MooseArray< Real > & JxW() const
Returns the reference to the transformed jacobian weights.
void reinitMortarElem(const Elem *elem)
reinitialize a mortar segment mesh element in order to get a proper JxW
Real _current_side_volume
Volume of the current side element.
std::vector< std::vector< Real > > _cached_residual_values
Values cached by calling cacheResidual() (the first vector is for TIME vs NONTIME) ...
virtual bool usesSecondPhiNeighbor() const =0
Whether or not this variable is actually using the shape function second derivatives.
virtual const FieldVariablePhiSecond & secondPhiNeighbor() const =0
Return the rank-2 tensor of second derivatives of the variable's shape functions on a neighboring ele...
void buildVectorFE(FEType type) const
Build Vector FEs with a type.
virtual const FieldVariablePhiGradient & gradPhiFaceNeighbor() const =0
Return the gradients of the variable's shape functions on a neighboring element face.
MooseArray< DualReal > _ad_curvatures
void reinitFVFace(const FaceInfo &fi)
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_lower
Vector FE objects for lower dimensional elements.
const MooseArray< Real > & JxWNeighbor() const
Returns the reference to the transformed jacobian weights on a current face.
std::vector< VectorValue< DualReal > > _ad_d2xyzdxideta_map
const Node * _current_node
The current node we are working with.
bool _building_helpers
Whether we are currently building the FE classes for the helpers.
virtual unsigned int numMatrixTags() const
The total number of tags.
OutputTools< Real >::VariablePhiGradient VariablePhiGradient
virtual MooseVariableScalar & getScalarVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a scalar variable with specified number.
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_face_neighbor
void addJacobian(GlobalDataKey)
Adds all local Jacobian to the global Jacobian matrices.
MooseArray< Point > _current_normals
The current Normal vectors at the quadrature points.
virtual const FieldVariablePhiValue & phiFace() const =0
Return the variable's shape functions on an element face.
void buildVectorFaceFE(FEType type) const
Build Vector FEs for a face with a type.
const SubdomainID ANY_BLOCK_ID
const VariablePhiGradient & gradPhiNeighbor(const MooseVariableField< Real > &) const
virtual const FieldVariablePhiGradient & gradPhi() const =0
Return the gradients of the variable's elemental shape functions.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void hasScalingVector()
signals this object that a vector containing variable scaling factors should be used when doing resid...
Generic class for solving transient nonlinear problems.
void setNeighborQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for neighbor integration.
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_dual_lower
void zeroCachedJacobian(GlobalDataKey)
Zero out previously-cached Jacobian rows.
std::vector< VectorValue< DualReal > > _ad_dxyzdzeta_map
void setCoordinateTransformation(const QBase *qrule, const Points &q_points, Coords &coord, SubdomainID sub_id)
void reinitNeighborFaceRef(const Elem *neighbor_elem, unsigned int neighbor_side, Real tolerance, const std::vector< Point > *const pts, const std::vector< Real > *const weights=nullptr)
Reinitialize FE data for the given neighbor_element on the given side with a given set of reference p...
bool usesGradPhiNeighbor() const
Whether or not this variable is actually using the shape function gradient.
void prepareVariableNonlocal(MooseVariableFieldBase *var)
void release()
Manually deallocates the data pointer.
std::map< FEType, std::unique_ptr< FEShapeData > > _fe_shape_data_face
virtual SparseMatrix< Number > & getMatrix(TagID tag)
Get a raw SparseMatrix.
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data
Shape function values, gradients, second derivatives for each vector FE type.
void addResidualScalar(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add residuals of all scalar variables for a set of tags onto the global residual vectors associated w...
MooseArray< Real > _coord_msm
The coordinate transformation coefficients evaluated on the quadrature points of the mortar segment m...
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Keg
bool _need_lower_d_elem_volume
Whether we need to compute the lower dimensional element volume.
bool _user_added_fe_face_of_helper_type
const std::vector< dof_id_type > & allDofIndices() const
Get all global dofindices for the variable.
void buildVectorNeighborFE(FEType type) const
Build Vector FEs for a neighbor with a type.
std::unique_ptr< ArbitraryQuadrature > arbitrary_vol
volume/elem (meshdim) custom points quadrature rule
void reinitFENeighbor(const Elem *neighbor, const std::vector< Point > &reference_points)
void initNonlocalCoupling()
Create pair of variables requiring nonlocal jacobian contributions.
std::map< unsigned int, FEBase * > _holder_fe_face_helper
Each dimension's helper objects.
void setLowerQRule(QBase *qrule, unsigned int dim)
Set the qrule to be used for lower dimensional integration.
const Elem *const & neighbor() const
Return the neighbor element.
std::map< unsigned int, std::map< FEType, FEVectorBase * > > _vector_fe_neighbor
void cacheJacobianNonlocal(GlobalDataKey)
Takes the values that are currently in _sub_Keg and appends them to the cached values.
Class for scalar variables (they are different).
IntRange< T > make_range(T beg, T end)
std::vector< std::unique_ptr< FEBase > > _unique_fe_face_neighbor_helper
void resize(unsigned int size)
Change the number of elements the array can store.
void computeSinglePointMapAD(const Elem *elem, const std::vector< Real > &qw, unsigned p, FEBase *fe)
compute the finite element reference-physical mapping quantities (such as JxW) with possible dependen...
void helpersRequestData()
request phi, dphi, xyz, JxW, etc.
std::vector< std::pair< MooseVariableFieldBase *, MooseVariableFieldBase * > > _cm_nonlocal_entry
Entries in the coupling matrix for field variables for nonlocal calculations.
void addJacobianOffDiagScalar(unsigned int ivar, GlobalDataKey)
Add Jacobians for a scalar variables with all other field variables into the global Jacobian matrices...
std::vector< DualReal > _ad_detadx_map
void buildLowerDFE(FEType type) const
Build FEs for a lower dimensional element with a type.
virtual unsigned int size() const override final
unsigned int _current_side
The current side of the selected element (valid only when working with sides)
void setCachedJacobian(GlobalDataKey)
Sets previously-cached Jacobian values via SparseMatrix::set() calls.
VariablePhiSecond _second_phi
std::map< unsigned int, FEBase * > _holder_fe_helper
Each dimension's helper objects.
void prepareVariable(MooseVariableFieldBase *var)
Used for preparing the dense residual and jacobian blocks for one particular variable.
const Elem * _current_neighbor_side_elem
The current side element of the ncurrent neighbor element.
void reinitElemFaceRef(const Elem *elem, unsigned int elem_side, Real tolerance, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)
Reinitialize FE data for the given element on the given side, optionally with a given set of referenc...
void copyFaceShapes(MooseVariableField< T > &v)
std::vector< dof_id_type > _extra_elem_ids
Extra element IDs.
void derivInsert(NumberArray< N, Real > &derivs, dof_id_type index, Real value)
Moose::CoordinateSystemType getCoordSystem(SubdomainID sid) const
void clearCachedQRules()
Set the cached quadrature rules to nullptr.
void addResidualNeighbor(GlobalDataKey, const std::vector< VectorTag > &vector_tags)
Add local neighbor residuals of all field variables for a set of tags onto the global residual vector...
bool _current_side_volume_computed
Boolean to indicate whether current element side volumes has been computed.
virtual bool isScalarVariable(unsigned int var_name) const
std::map< FEType, ADTemplateVariablePhiGradient< RealVectorValue > > _ad_vector_grad_phi_data_face
std::vector< VectorValue< DualReal > > _ad_d2xyzdeta2_map
Moose::VectorTagType _type
The type of the vector tag.
void jacobianBlockNeighborUsed(TagID tag, unsigned int ivar, unsigned int jvar, bool used)
Sets whether or not neighbor Jacobian coupling between ivar and jvar is used to the value used...
std::unique_ptr< ArbitraryQuadrature > arbitrary_face
area/face (meshdim-1) custom points quadrature rule
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
virtual const std::vector< dof_id_type > & dofIndicesLower() const =0
Get dof indices for the current lower dimensional element (this is meaningful when performing mortar ...
void prepareBlock(unsigned int ivar, unsigned jvar, const std::vector< dof_id_type > &dof_indices)
std::map< unsigned int, std::map< FEType, FEBase * > > _fe_face_neighbor
void addJacobianNeighbor(GlobalDataKey)
Add ElementNeighbor, NeighborElement, and NeighborNeighbor portions of the Jacobian for compute objec...
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_lower_used
Flag that indicates if the jacobian block for the lower dimensional element was used.
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
VectorVariablePhiDivergence _div_phi
void cacheJacobianBlock(DenseMatrix< Number > &jac_block, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, Real scaling_factor, LocalDataKey, TagID tag)
Cache a local Jacobian block with the provided rows (idof_indices) and columns (jdof_indices) for eve...
MooseVariableFieldBase & getVariable(THREAD_ID tid, const std::string &var_name) const
Gets a reference to a variable of with specified name.
bool _custom_mortar_qrule
Flag specifying whether a custom quadrature rule has been specified for mortar segment mesh...
const unsigned int & side() const
Returns the current side.
void computeCurrentElemVolume()
void cacheJacobianBlockNonzero(DenseMatrix< Number > &jac_block, const MooseVariableBase &ivar, const MooseVariableBase &jvar, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices, TagID tag)
Push non-zeros of a local Jacobian block with proper scaling into cache for a certain tag...
Real _current_neighbor_lower_d_elem_volume
The current neighboring lower dimensional element volume.
bool _need_neighbor_lower_d_elem_volume
Whether we need to compute the neighboring lower dimensional element volume.
MooseArray< Point > _current_q_points
The current list of quadrature points.
bool _user_added_fe_of_helper_type
Whether user code requested a FEType the same as our _helper_type.
Moose::CoordinateSystemType _coord_type
The coordinate system.
void cacheJacobianMortar(GlobalDataKey)
Cache all portions of the Jacobian, e.g.
std::unique_ptr< ArbitraryQuadrature > neighbor
area/face (meshdim-1) custom points quadrature rule for DG
void addJacobianBlock(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, std::vector< dof_id_type > &dof_indices, GlobalDataKey, TagID tag)
Adds element matrix for ivar rows and jvar columns to the global Jacobian matrix. ...
void buildVectorFaceNeighborFE(FEType type) const
Build Vector FEs for a neighbor face with a type.
void buildVectorDualLowerDFE(FEType type) const
void bumpVolumeQRuleOrder(Order volume_order, SubdomainID block)
Increases the element/volume quadrature order for the specified mesh block if and only if the current...
Storage for all of the information pretaining to a vector tag.
const Node *const & node() const
Returns the reference to the node.
std::vector< std::vector< std::vector< unsigned char > > > _jacobian_block_neighbor_used
Flag that indicates if the jacobian block for neighbor was used.
std::map< unsigned int, FEBase * > _holder_fe_lower_helper
helper object for transforming coordinates for lower dimensional element quadrature points ...
void computeADFace(const Elem &elem, const unsigned int side)
compute AD things on an element face
std::map< FEType, std::unique_ptr< VectorFEShapeData > > _vector_fe_shape_data_neighbor
std::vector< DualReal > _ad_detady_map
void setResidual(NumericVector< Number > &residual, GlobalDataKey, const VectorTag &vector_tag)
Sets local residuals of all field variables to the global residual vector for a tag.
std::vector< std::vector< dof_id_type > > _cached_jacobian_cols
Column where the corresponding cached value should go.
void addCachedJacobian(GlobalDataKey)
Adds the values that have been cached by calling cacheJacobian() and or cacheJacobianNeighbor() to th...
void addJacobianBlockTags(SparseMatrix< Number > &jacobian, unsigned int ivar, unsigned int jvar, const DofMap &dof_map, std::vector< dof_id_type > &dof_indices, GlobalDataKey, const std::set< TagID > &tags)
Add element matrix for ivar rows and jvar columns to the global Jacobian matrix for given tags...
DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
MooseArray< Point > _current_q_points_face
The current quadrature points on a face.
void cacheResidualBlock(std::vector< Real > &cached_residual_values, std::vector< dof_id_type > &cached_residual_rows, DenseVector< Number > &res_block, const std::vector< dof_id_type > &dof_indices, const std::vector< Real > &scaling_factor, bool is_nodal)
Push a local residual block with proper scaling into cache.
void prepareBlockNonlocal(unsigned int ivar, unsigned jvar, const std::vector< dof_id_type > &idof_indices, const std::vector< dof_id_type > &jdof_indices)
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector by name.
void bumpAllQRuleOrder(Order order, SubdomainID block)
Increases the element/volume and face/area quadrature orders for the specified mesh block if and only...
const VariablePhiSecond & secondPhiFaceNeighbor(const MooseVariableField< Real > &) const
const Elem * _current_side_elem
The current "element" making up the side we are currently on.
bool _block_diagonal_matrix
Will be true if our preconditioning matrix is a block-diagonal matrix. Which means that we can take s...
MooseArray< DualReal > _ad_coord
The AD version of the current coordinate transformation coefficients.
bool _need_JxW_neighbor
Flag to indicate that JxW_neighbor is needed.
std::map< FEType, bool > _need_curl
virtual const VectorTag & getVectorTag(const TagID tag_id) const
Get a VectorTag from a TagID.
void copyNeighborShapes(MooseVariableField< T > &v)
void cacheJacobianCoupledVarPair(const MooseVariableBase &ivar, const MooseVariableBase &jvar)
Caches element matrix for ivar rows and jvar columns.
MooseArray< Point > _current_q_points_face_neighbor
The current quadrature points on the neighbor face.
VectorVariablePhiCurl _vector_curl_phi_face
Key structure for APIs adding/caching local element residuals/Jacobians.
SubdomainID _current_subdomain_id
The current subdomain ID.
virtual numeric_index_type n() const =0
ElemSideBuilder _current_side_elem_builder
In place side element builder for _current_side_elem.
QBase * _current_qrule_face
quadrature rule used on faces
std::map< FEType, FEVectorBase * > _current_vector_fe_neighbor
The "neighbor" vector fe object that matches the current elem.
void clearCachedResiduals(GlobalDataKey)
Clears all of the residuals in _cached_residual_rows and _cached_residual_values. ...
virtual const FieldVariablePhiGradient & gradPhiFace() const =0
Return the gradients of the variable's shape functions on an element face.
const Node * _current_neighbor_node
The current neighboring node we are working with.
void clearCachedJacobian()
Clear any currently cached jacobians.
const VariablePhiGradient & gradPhiFace() const
std::vector< std::vector< std::vector< DenseMatrix< Number > > > > _sub_Knl
dprimary/dlower (or dneighbor/dlower)
std::vector< VectorValue< DualReal > > _ad_dxyzdxi_map
AD quantities.
MooseVariableField< T > & getActualFieldVariable(THREAD_ID tid, const std::string &var_name)
Returns a field variable pointer - this includes finite volume variables.
void reinitFEFace(const Elem *elem, unsigned int side)
Just an internal helper function to reinit the face FE objects.
void computeFaceMap(const Elem &elem, const unsigned int side, const std::vector< Real > &qw)
void reinitDual(const Elem *elem, const std::vector< Point > &pts, const std::vector< Real > &JxW)
Reintialize dual basis coefficients based on a customized quadrature rule.
Key structure for APIs manipulating global vectors/matrices.