https://mooseframework.inl.gov
EqualValueEmbeddedConstraint.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 // MOOSE includes
11 #include "EigenADReal.h"
13 #include "FEProblem.h"
14 #include "DisplacedProblem.h"
15 #include "AuxiliarySystem.h"
16 #include "SystemBase.h"
17 #include "Assembly.h"
18 #include "MooseMesh.h"
19 #include "AddVariableAction.h"
20 #include "MooseEnum.h"
21 
22 #include "libmesh/sparse_matrix.h"
23 
26 
27 template <bool is_ad>
30 {
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");
37  params.addParam<MooseEnum>(
38  "formulation", formulation, "Formulation used to enforce the constraint");
39  params.addRequiredParam<Real>(
40  "penalty",
41  "Penalty parameter used in constraint enforcement for kinematic and penalty formulations.");
42 
43  return params;
44 }
45 
46 template <bool is_ad>
48  const InputParameters & parameters)
49  : GenericNodeElemConstraint<is_ad>(parameters),
50  _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
51  _fe_problem(*parameters.get<FEProblem *>("_fe_problem")),
52  _formulation(this->template getParam<MooseEnum>("formulation").template getEnum<Formulation>()),
53  _penalty(this->template getParam<Real>("penalty")),
54  _residual_copy(_sys.residualGhosted())
55 {
58  if constexpr (is_ad)
59  {
61  this->paramError("formulation", "AD constraints cannot be used with KINEMATIC formulation.");
62  }
63 }
64 
65 template <bool is_ad>
66 void
68 {
69  // get mesh pointLocator
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};
73 
74  // secondary id and primary id
75  dof_id_type sid, mid;
76 
77  // prepare _secondary_to_primary_map
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)))
82  {
83  for (auto & sn : elem->node_ref_range())
84  {
85  sid = sn.id();
86  if (_secondary_to_primary_map.find(sid) == _secondary_to_primary_map.end())
87  {
88  // primary element
89  const Elem * me = pointLocator->operator()(sn, &allowed_subdomains);
90  if (me != NULL)
91  {
92  mid = me->id();
93  _secondary_to_primary_map.insert(std::pair<dof_id_type, dof_id_type>(sid, mid));
94  _subproblem.addGhostedElem(mid);
95  }
96  }
97  }
98  }
99 }
100 
101 template <bool is_ad>
102 bool
104 {
105  // primary element
106  auto it = _secondary_to_primary_map.find(_current_node->id());
107 
108  if (it != _secondary_to_primary_map.end())
109  {
110  const Elem * primary_elem = _mesh.elemPtr(it->second);
111  std::vector<Point> points = {*_current_node};
112 
113  // reinit variables on the primary element at the secondary point
114  _fe_problem.setNeighborSubdomainID(primary_elem, 0);
115  _fe_problem.reinitNeighborPhys(primary_elem, points, 0);
116 
117  reinitConstraint();
118 
119  return true;
120  }
121  return false;
122 }
123 
124 template <bool is_ad>
125 void
127 {
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);
131 
132  switch (_formulation)
133  {
134  case Formulation::KINEMATIC:
135  _constraint_residual = -_residual_copy(dof_number);
136  break;
137 
138  case Formulation::PENALTY:
139  _constraint_residual = _penalty * (_u_secondary[0] - _u_primary[0]);
140  break;
141 
142  default:
143  mooseError("Invalid formulation");
144  break;
145  }
146 }
147 
148 template <bool is_ad>
149 Real
151 {
152  return MetaPhysicL::raw_value(_u_secondary[_qp]);
153 }
154 
155 template <bool is_ad>
158 {
159  GenericReal<is_ad> resid = _constraint_residual;
160 
161  switch (type)
162  {
163  case Moose::Secondary:
164  {
165  if (_formulation == Formulation::KINEMATIC)
166  {
167  GenericReal<is_ad> pen_force = _penalty * (_u_secondary[_qp] - _u_primary[_qp]);
168  resid += pen_force;
169  }
170  return _test_secondary[_i][_qp] * resid;
171  }
172 
173  case Moose::Primary:
174  return _test_primary[_i][_qp] * -resid;
175  }
176 
177  return 0.0;
178 }
179 
180 template <bool is_ad>
181 Real
183 {
184  mooseAssert(!is_ad,
185  "In ADEqualValueEmbeddedConstraint, computeQpJacobian should not be called. "
186  "Check computeJacobian implementation.");
187 
188  unsigned int sys_num = _sys.number();
189  const Real penalty = MetaPhysicL::raw_value(_penalty);
190  Real curr_jac, secondary_jac;
191 
192  switch (type)
193  {
195  switch (_formulation)
196  {
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];
203  default:
204  mooseError("Invalid formulation");
205  }
206 
208  switch (_formulation)
209  {
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];
214  default:
215  mooseError("Invalid formulation");
216  }
217 
219  switch (_formulation)
220  {
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];
227  default:
228  mooseError("Invalid formulation");
229  }
230 
232  switch (_formulation)
233  {
234  case Formulation::KINEMATIC:
235  return 0.0;
236  case Formulation::PENALTY:
237  return _test_primary[_i][_qp] * penalty * _phi_primary[_j][_qp];
238  default:
239  mooseError("Invalid formulation");
240  }
241 
242  default:
243  mooseError("Unsupported type");
244  break;
245  }
246  return 0.0;
247 }
248 
249 template <bool is_ad>
250 Real
252  Moose::ConstraintJacobianType type, unsigned int /*jvar*/)
253 {
254  mooseAssert(!is_ad,
255  "In ADEqualValueEmbeddedConstraint, computeQpOffDiagJacobian should not be called. "
256  "Check computeJacobian implementation.");
257 
258  Real curr_jac, secondary_jac;
259  unsigned int sys_num = _sys.number();
260 
261  switch (type)
262  {
264  curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
265  _connected_dof_indices[_j]);
266  return -curr_jac;
267 
269  return 0.0;
270 
272  switch (_formulation)
273  {
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:
279  return 0.0;
280  default:
281  mooseError("Invalid formulation");
282  }
283 
285  return 0.0;
286 
287  default:
288  mooseError("Unsupported type");
289  break;
290  }
291 
292  return 0.0;
293 }
294 
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
Definition: MooseTypes.h:648
ConstraintType
Definition: MooseTypes.h:757
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.
Definition: FEProblem.h:20
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
T * get(const std::unique_ptr< T > &u)
The MooseUtils::get() specializations are used to support making forwards-compatible code changes fro...
Definition: MooseUtils.h:1155
virtual void prepareSecondaryToPrimaryMap() override
prepare the _secondary_to_primary_map (see NodeElemConstraintBase)
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
bool _overwrite_secondary_residual
Whether or not the secondary&#39;s residual should be overwritten.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
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...
Definition: MooseEnum.h:33
void paramError(const std::string &param, 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
ConstraintJacobianType
Definition: MooseTypes.h:796
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an optional parameter and a documentation string to the InputParameters object...
EqualValueEmbeddedConstraintTempl(const InputParameters &parameters)
static InputParameters validParams()
uint8_t dof_id_type