22 #include "libmesh/sparse_matrix.h" 33 params.
addClassDescription(
"This is a constraint enforcing overlapping portions of two blocks to " 34 "have the same variable value");
35 params.
set<
bool>(
"use_displaced_mesh") =
false;
36 MooseEnum formulation(
"kinematic penalty",
"kinematic");
38 "formulation", formulation,
"Formulation used to enforce the constraint");
41 "Penalty parameter used in constraint enforcement for kinematic and penalty formulations.");
50 _displaced_problem(parameters.
get<
FEProblemBase *>(
"_fe_problem_base")->getDisplacedProblem()),
51 _fe_problem(*parameters.
get<
FEProblem *>(
"_fe_problem")),
53 _penalty(this->template getParam<
Real>(
"penalty")),
54 _residual_copy(_sys.residualGhosted())
61 this->
paramError(
"formulation",
"AD constraints cannot be used with KINEMATIC formulation.");
70 std::unique_ptr<libMesh::PointLocatorBase> pointLocator = _mesh.getPointLocator();
71 pointLocator->enable_out_of_mesh_mode();
72 const std::set<subdomain_id_type> allowed_subdomains{_primary};
78 std::set<dof_id_type> unique_secondary_node_ids;
79 const MeshBase & meshhelper = _mesh.getMesh();
80 for (
const auto & elem :
as_range(meshhelper.active_subdomain_elements_begin(_secondary),
81 meshhelper.active_subdomain_elements_end(_secondary)))
83 for (
auto & sn : elem->node_ref_range())
86 if (_secondary_to_primary_map.find(sid) == _secondary_to_primary_map.end())
89 const Elem * me = pointLocator->operator()(sn, &allowed_subdomains);
93 _secondary_to_primary_map.insert(std::pair<dof_id_type, dof_id_type>(sid, mid));
94 _subproblem.addGhostedElem(mid);
101 template <
bool is_ad>
106 auto it = _secondary_to_primary_map.find(_current_node->id());
108 if (it != _secondary_to_primary_map.end())
110 const Elem * primary_elem = _mesh.elemPtr(it->second);
111 std::vector<Point> points = {*_current_node};
114 _fe_problem.setNeighborSubdomainID(primary_elem, 0);
115 _fe_problem.reinitNeighborPhys(primary_elem, points, 0);
124 template <
bool is_ad>
128 const Node * node = _current_node;
129 unsigned int sys_num = _sys.number();
130 dof_id_type dof_number = node->dof_number(sys_num, _var.number(), 0);
132 switch (_formulation)
134 case Formulation::KINEMATIC:
135 _constraint_residual = -_residual_copy(dof_number);
138 case Formulation::PENALTY:
139 _constraint_residual = _penalty * (_u_secondary[0] - _u_primary[0]);
148 template <
bool is_ad>
155 template <
bool is_ad>
165 if (_formulation == Formulation::KINEMATIC)
170 return _test_secondary[_i][_qp] * resid;
174 return _test_primary[_i][_qp] * -resid;
180 template <
bool is_ad>
185 "In ADEqualValueEmbeddedConstraint, computeQpJacobian should not be called. " 186 "Check computeJacobian implementation.");
188 unsigned int sys_num = _sys.number();
190 Real curr_jac, secondary_jac;
195 switch (_formulation)
197 case Formulation::KINEMATIC:
198 curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
199 _connected_dof_indices[_j]);
200 return -curr_jac + _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
201 case Formulation::PENALTY:
202 return _phi_secondary[_j][_qp] * penalty * _test_secondary[_i][_qp];
208 switch (_formulation)
210 case Formulation::KINEMATIC:
211 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
212 case Formulation::PENALTY:
213 return -_phi_primary[_j][_qp] * penalty * _test_secondary[_i][_qp];
219 switch (_formulation)
221 case Formulation::KINEMATIC:
222 secondary_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
223 _connected_dof_indices[_j]);
224 return secondary_jac * _test_primary[_i][_qp];
225 case Formulation::PENALTY:
226 return -_phi_secondary[_j][_qp] * penalty * _test_primary[_i][_qp];
232 switch (_formulation)
234 case Formulation::KINEMATIC:
236 case Formulation::PENALTY:
237 return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
249 template <
bool is_ad>
255 "In ADEqualValueEmbeddedConstraint, computeQpOffDiagJacobian should not be called. " 256 "Check computeJacobian implementation.");
258 Real curr_jac, secondary_jac;
259 unsigned int sys_num = _sys.number();
264 curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
265 _connected_dof_indices[_j]);
272 switch (_formulation)
274 case Formulation::KINEMATIC:
275 secondary_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
276 _connected_dof_indices[_j]);
277 return secondary_jac * _test_primary[_i][_qp];
278 case Formulation::PENALTY:
static InputParameters validParams()
virtual GenericReal< is_ad > computeQpResidual(Moose::ConstraintType type) override
This is the virtual method that derived classes should override for computing the residual...
Moose::GenericType< Real, is_ad > GenericReal
virtual Real computeQpJacobian(Moose::ConstraintJacobianType type) override
This is the virtual method that derived classes should override for computing the Jacobian...
A EqualValueEmbeddedConstraint forces the value of a variable to be the same on overlapping portion o...
virtual void reinitConstraint()
Prepare the residual contribution of the current constraint required to enforce it based on the speci...
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
T * get(const std::unique_ptr< T > &u)
The MooseUtils::get() specializations are used to support making forwards-compatible code changes fro...
virtual void prepareSecondaryToPrimaryMap() override
prepare the _secondary_to_primary_map (see NodeElemConstraintBase)
bool _overwrite_secondary_residual
Whether or not the secondary's residual should be overwritten.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
bool shouldApply() override
Whether or not this constraint should be applied.
static MooseEnum getNonlinearVariableOrders()
Get the possible variable orders.
registerMooseObject("MooseApp", EqualValueEmbeddedConstraint)
enum EqualValueEmbeddedConstraintTempl::Formulation _formulation
Formulation
Formulations, currently only supports KINEMATIC and PENALTY.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
void paramError(const std::string ¶m, Args... args) const
Emits an error prefixed with the file and line number of the given param (from the input file) along ...
virtual Real computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, unsigned int jvar) override
This is the virtual method that derived classes should override for computing the off-diag Jacobian...
virtual Real computeQpSecondaryValue() override
Compute the value the secondary node should have at the beginning of a timestep.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
EqualValueEmbeddedConstraintTempl(const InputParameters ¶meters)
static InputParameters validParams()