14 #include "metaphysicl/dualsemidynamicsparsenumberarray.h" 15 #include "metaphysicl/parallel_dualnumber.h" 16 #include "metaphysicl/parallel_dynamic_std_array_wrapper.h" 17 #include "metaphysicl/parallel_semidynamicsparsenumberarray.h" 28 const auto & disp_x_name = ret.
get<std::vector<VariableName>>(
"disp_x");
29 if (disp_x_name.size() != 1)
30 mooseError(
"We require that the disp_x parameter have exactly one coupled name");
33 ret.
set<VariableName>(
"secondary_variable") = disp_x_name[0];
34 ret.
set<VariableName>(
"primary_variable") = disp_x_name[0];
45 params.
addParam<
Real>(
"c_t", 1e0,
"Numerical parameter for tangential constraints");
49 "Minimum value of contact pressure that will trigger frictional enforcement");
51 "mu",
"mu > 0",
"The friction coefficient for the Coulomb friction law");
58 _c_t(getParam<
Real>(
"c_t")),
59 _secondary_x_dot(adCoupledDot(
"disp_x")),
60 _primary_x_dot(adCoupledNeighborValueDot(
"disp_x")),
61 _secondary_y_dot(adCoupledDot(
"disp_y")),
62 _primary_y_dot(adCoupledNeighborValueDot(
"disp_y")),
63 _secondary_z_dot(_has_disp_z ? &adCoupledDot(
"disp_z") : nullptr),
64 _primary_z_dot(_has_disp_z ? &adCoupledNeighborValueDot(
"disp_z") : nullptr),
65 _mu(getParam<
Real>(
"mu")),
66 _epsilon(getParam<
Real>(
"epsilon"))
78 const auto & secondary_ip_lowerd_map =
82 std::array<ADReal, 3> primary_disp_dot{
84 std::array<ADReal, 3> secondary_disp_dot{
90 const ADReal & prim_x_dot = primary_disp_dot[0];
91 const ADReal & prim_y_dot = primary_disp_dot[1];
92 const ADReal * prim_z_dot =
nullptr;
94 prim_z_dot = &primary_disp_dot[2];
96 const ADReal & sec_x_dot = secondary_disp_dot[0];
97 const ADReal & sec_y_dot = secondary_disp_dot[1];
98 const ADReal * sec_z_dot =
nullptr;
100 sec_z_dot = &secondary_disp_dot[2];
106 relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, *sec_z_dot - *prim_z_dot};
108 relative_velocity = {sec_x_dot - prim_x_dot, sec_y_dot - prim_y_dot, 0.0};
120 const DofObject *
const dof =
151 const DofObject *
const dof = pr.first;
170 const std::unordered_set<const Node *> & inactive_lm_nodes)
180 const DofObject *
const dof = pr.first;
183 if ((inactive_lm_nodes.find(static_cast<const Node *>(dof)) != inactive_lm_nodes.end()) ||
200 const DofObject *
const dof)
208 const Real scaling_factor_x =
_lm_vars[0]->scalingFactor();
209 const Real scaling_factor_y =
_lm_vars[1]->scalingFactor();
210 Real scaling_factor_z = 1;
225 scaling_factor_z =
_lm_vars[2]->scalingFactor();
228 ADReal normal_pressure_value =
230 ADReal tangential_pressure_value =
233 ADReal tangential_pressure_value_dir;
244 ADReal normal_dof_residual = std::min(normal_pressure_value, weighted_gap *
c);
245 ADReal tangential_dof_residual;
246 ADReal tangential_dof_residual_dir;
253 if (normal_pressure_value.value() <
_epsilon)
254 tangential_dof_residual = tangential_pressure_value;
258 std::max(
_mu * (normal_pressure_value +
c * weighted_gap),
259 std::abs(tangential_pressure_value + c_t * tangential_vel *
_dt)) *
260 tangential_pressure_value;
261 const auto term_2 =
_mu * std::max(0.0, normal_pressure_value +
c * weighted_gap) *
262 (tangential_pressure_value + c_t * tangential_vel *
_dt);
264 tangential_dof_residual = term_1 - term_2;
269 if (normal_pressure_value.value() <
_epsilon)
271 tangential_dof_residual = tangential_pressure_value;
272 tangential_dof_residual_dir = tangential_pressure_value_dir;
276 const Real epsilon_sqrt = 1.0e-48;
278 const auto lamdba_plus_cg = normal_pressure_value +
c * weighted_gap;
280 std::array<ADReal, 2> lambda_t_plus_ctu;
284 const auto term_1_x =
285 std::max(
_mu * lamdba_plus_cg,
286 std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
287 lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
288 tangential_pressure_value;
290 const auto term_1_y =
291 std::max(
_mu * lamdba_plus_cg,
292 std::sqrt(lambda_t_plus_ctu[0] * lambda_t_plus_ctu[0] +
293 lambda_t_plus_ctu[1] * lambda_t_plus_ctu[1] + epsilon_sqrt)) *
294 tangential_pressure_value_dir;
296 const auto term_2_x =
_mu * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[0];
298 const auto term_2_y =
_mu * std::max(0.0, lamdba_plus_cg) * lambda_t_plus_ctu[1];
300 tangential_dof_residual = term_1_x - term_2_x;
301 tangential_dof_residual_dir = term_1_y - term_2_y;
324 unsigned int component_normal = 0;
327 const Real threshold_for_Jacobian =
_has_disp_z ? 1.0 / std::sqrt(3.0) : 1.0 / std::sqrt(2.0);
329 if (std::abs(ny) > threshold_for_Jacobian)
330 component_normal = 1;
331 else if (std::abs(nz) > threshold_for_Jacobian)
332 component_normal = 2;
336 std::array<ADReal, 1>{{normal_dof_residual}},
337 std::array<dof_id_type, 1>{{component_normal == 0
339 : (component_normal == 1 ? dof_index_y : dof_index_z)}},
340 component_normal == 0 ? scaling_factor_x
341 : (component_normal == 1 ? scaling_factor_y : scaling_factor_z));
345 std::array<ADReal, 1>{{tangential_dof_residual}},
346 std::array<dof_id_type, 1>{
347 {(component_normal == 0 || component_normal == 2) ? dof_index_y : dof_index_x}},
348 (component_normal == 0 || component_normal == 2) ? scaling_factor_y : scaling_factor_x);
353 std::array<ADReal, 1>{{tangential_dof_residual_dir}},
354 std::array<dof_id_type, 1>{
355 {(component_normal == 0 || component_normal == 1) ? dof_index_z : dof_index_x}},
356 (component_normal == 0 || component_normal == 1) ? scaling_factor_z : scaling_factor_x);
virtual const NumericVector< Number > *const & currentSolution() const=0
void addResidualsAndJacobian(Assembly &assembly, const Residuals &residuals, const Indices &dof_indices, Real scaling_factor)
void mooseError(Args &&... args)
std::map< unsigned int, unsigned int > getPrimaryIpToLowerElementMap(const Elem &primary_elem, const Elem &primary_elem_ip, const Elem &lower_secondary_elem) const
Elem const *const & _lower_primary_elem
const std::vector< Real > & _JxW_msm
const Parallel::Communicator & _communicator
DualNumber< Real, DNDerivativeType, true > ADReal
const VariableTestValue & _test
Elem const *const & _lower_secondary_elem
const AutomaticMortarGeneration & amg() const
static void trimInteriorNodeDerivatives(const std::map< unsigned int, unsigned int > &primary_ip_lowerd_map, const Variables &moose_var, DualNumbers &ad_vars, const bool is_secondary)
std::array< MooseUtils::SemidynamicVector< Point, 9 >, 2 > getNodalTangents(const Elem &secondary_elem) const
std::map< unsigned int, unsigned int > getSecondaryIpToLowerElementMap(const Elem &lower_secondary_elem) const
unsigned int number() const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void derivInsert(SemiDynamicSparseNumberArray< Real, libMesh::dof_id_type, NWrapper< N >> &derivs, libMesh::dof_id_type index, Real value)
const MooseArray< Real > & _coord
processor_id_type processor_id() const