www.mooseframework.org
EqualValueEmbeddedConstraint.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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
12 #include "FEProblem.h"
13 #include "DisplacedProblem.h"
14 #include "AuxiliarySystem.h"
15 #include "SystemBase.h"
16 #include "Assembly.h"
17 #include "MooseMesh.h"
18 #include "Executioner.h"
19 #include "AddVariableAction.h"
20 
21 #include "libmesh/string_to_enum.h"
22 #include "libmesh/sparse_matrix.h"
23 
25 
26 template <>
29 {
32  params.addClassDescription("This is a constraint enforcing overlapping portions of two blocks to "
33  "have the same variable value");
34  params.set<bool>("use_displaced_mesh") = false;
35  MooseEnum formulation("kinematic penalty", "kinematic");
36  params.addParam<MooseEnum>(
37  "formulation", formulation, "Formulation used to enforce the constraint");
38  params.addRequiredParam<Real>(
39  "penalty",
40  "Penalty parameter used in constraint enforcement for kinematic and penalty formulations.");
41 
42  return params;
43 }
44 
46  : NodeElemConstraint(parameters),
47  _displaced_problem(parameters.get<FEProblemBase *>("_fe_problem_base")->getDisplacedProblem()),
48  _fe_problem(*parameters.get<FEProblem *>("_fe_problem")),
49  _formulation(getParam<MooseEnum>("formulation").getEnum<Formulation>()),
50  _penalty(getParam<Real>("penalty")),
51  _residual_copy(_sys.residualGhosted())
52 {
55 }
56 
57 void
59 {
60  // get mesh pointLocator
61  std::unique_ptr<PointLocatorBase> pointLocator = _mesh.getPointLocator();
62  pointLocator->enable_out_of_mesh_mode();
63  const std::set<subdomain_id_type> allowed_subdomains{_master};
64 
65  // slave id and master id
66  dof_id_type sid, mid;
67 
68  // prepare _slave_to_master_map
69  std::set<dof_id_type> unique_slave_node_ids;
70  const MeshBase & meshhelper = _mesh.getMesh();
71  for (const auto & elem : as_range(meshhelper.active_subdomain_elements_begin(_slave),
72  meshhelper.active_subdomain_elements_end(_slave)))
73  {
74  for (auto & sn : elem->node_ref_range())
75  {
76  sid = sn.id();
77  if (_slave_to_master_map.find(sid) == _slave_to_master_map.end())
78  {
79  // master element
80  const Elem * me = pointLocator->operator()(sn, &allowed_subdomains);
81  if (me != NULL)
82  {
83  mid = me->id();
84  _slave_to_master_map.insert(std::pair<dof_id_type, dof_id_type>(sid, mid));
86  }
87  }
88  }
89  }
90 }
91 
92 bool
94 {
95  // master element
96  auto it = _slave_to_master_map.find(_current_node->id());
97 
98  if (it != _slave_to_master_map.end())
99  {
100  const Elem * master_elem = _mesh.elemPtr(it->second);
101  std::vector<Point> points = {*_current_node};
102 
103  // reinit variables on the master element at the slave point
104  _fe_problem.setNeighborSubdomainID(master_elem, 0);
105  _fe_problem.reinitNeighborPhys(master_elem, points, 0);
106 
108 
109  return true;
110  }
111  return false;
112 }
113 
114 void
116 {
117  const Node * node = _current_node;
118  unsigned int sys_num = _sys.number();
119  dof_id_type dof_number = node->dof_number(sys_num, _var.number(), 0);
120 
121  switch (_formulation)
122  {
124  _constraint_residual = -_residual_copy(dof_number);
125  break;
126 
129  break;
130 
131  default:
132  mooseError("Invalid formulation");
133  break;
134  }
135 }
136 
137 Real
139 {
140  return _u_slave[_qp];
141 }
142 
143 Real
145 {
146  Real resid = _constraint_residual;
147 
148  switch (type)
149  {
150  case Moose::Slave:
151  {
153  {
154  Real pen_force = _penalty * (_u_slave[_qp] - _u_master[_qp]);
155  resid += pen_force;
156  }
157  return _test_slave[_i][_qp] * resid;
158  }
159 
160  case Moose::Master:
161  return _test_master[_i][_qp] * -resid;
162  }
163 
164  return 0.0;
165 }
166 
167 Real
169 {
170  unsigned int sys_num = _sys.number();
171  const Real penalty = _penalty;
172  Real curr_jac, slave_jac;
173 
174  switch (type)
175  {
176  case Moose::SlaveSlave:
177  switch (_formulation)
178  {
180  curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
182  return -curr_jac + _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
184  return _phi_slave[_j][_qp] * penalty * _test_slave[_i][_qp];
185  default:
186  mooseError("Invalid formulation");
187  }
188 
189  case Moose::SlaveMaster:
190  switch (_formulation)
191  {
193  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
195  return -_phi_master[_j][_qp] * penalty * _test_slave[_i][_qp];
196  default:
197  mooseError("Invalid formulation");
198  }
199 
200  case Moose::MasterSlave:
201  switch (_formulation)
202  {
204  slave_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
206  return slave_jac * _test_master[_i][_qp];
208  return -_phi_slave[_j][_qp] * penalty * _test_master[_i][_qp];
209  default:
210  mooseError("Invalid formulation");
211  }
212 
213  case Moose::MasterMaster:
214  switch (_formulation)
215  {
217  return 0.0;
219  return _test_master[_i][_qp] * penalty * _phi_master[_j][_qp];
220  default:
221  mooseError("Invalid formulation");
222  }
223 
224  default:
225  mooseError("Unsupported type");
226  break;
227  }
228  return 0.0;
229 }
230 
231 Real
233  unsigned int /*jvar*/)
234 {
235  Real curr_jac, slave_jac;
236  unsigned int sys_num = _sys.number();
237 
238  switch (type)
239  {
240  case Moose::SlaveSlave:
241  curr_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
243  return -curr_jac;
244 
245  case Moose::SlaveMaster:
246  return 0.0;
247 
248  case Moose::MasterSlave:
249  switch (_formulation)
250  {
252  slave_jac = (*_jacobian)(_current_node->dof_number(sys_num, _var.number(), 0),
254  return slave_jac * _test_master[_i][_qp];
256  return 0.0;
257  default:
258  mooseError("Invalid formulation");
259  }
260 
261  case Moose::MasterMaster:
262  return 0.0;
263 
264  default:
265  mooseError("Unsupported type");
266  break;
267  }
268 
269  return 0.0;
270 }
271 
272 void
274 {
276 
277  DenseMatrix<Number> & Knn =
279 
280  _Kee.resize(_test_slave.size(), _connected_dof_indices.size());
281 
282  for (_i = 0; _i < _test_slave.size(); _i++)
283  // Loop over the connected dof indices so we can get all the jacobian contributions
284  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
286 
287  DenseMatrix<Number> & Ken =
289  if (Ken.m() && Ken.n())
290  for (_i = 0; _i < _test_slave.size(); _i++)
291  for (_j = 0; _j < _phi_master.size(); _j++)
293 
294  _Kne.resize(_test_master.size(), _connected_dof_indices.size());
295  for (_i = 0; _i < _test_master.size(); _i++)
296  // Loop over the connected dof indices so we can get all the jacobian contributions
297  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
299 
300  if (Knn.m() && Knn.n())
301  for (_i = 0; _i < _test_master.size(); _i++)
302  for (_j = 0; _j < _phi_master.size(); _j++)
304 }
305 
306 void
308 {
310 
311  _Kee.resize(_test_slave.size(), _connected_dof_indices.size());
312 
313  DenseMatrix<Number> & Knn =
315 
316  for (_i = 0; _i < _test_slave.size(); _i++)
317  // Loop over the connected dof indices so we can get all the jacobian contributions
318  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
320 
321  DenseMatrix<Number> & Ken =
323  for (_i = 0; _i < _test_slave.size(); _i++)
324  for (_j = 0; _j < _phi_master.size(); _j++)
326 
327  _Kne.resize(_test_master.size(), _connected_dof_indices.size());
328  if (_Kne.m() && _Kne.n())
329  for (_i = 0; _i < _test_master.size(); _i++)
330  // Loop over the connected dof indices so we can get all the jacobian contributions
331  for (_j = 0; _j < _connected_dof_indices.size(); _j++)
333 
334  for (_i = 0; _i < _test_master.size(); _i++)
335  for (_j = 0; _j < _phi_master.size(); _j++)
337 }
338 
339 void
341 {
343 
345 
346  dof_id_type current_node_var_dof_index = _sys.getVariable(0, var_num).nodalDofIndex();
347 
348  // Fill up _phi_slave so that it is 1 when j corresponds to the dof associated with this node
349  // and 0 for every other dof
350  // This corresponds to evaluating all of the connected shape functions at _this_ node
351  _qp = 0;
352  for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
353  {
354  _phi_slave[j].resize(1);
355 
356  if (_connected_dof_indices[j] == current_node_var_dof_index)
357  _phi_slave[j][_qp] = 1.0;
358  else
359  _phi_slave[j][_qp] = 0.0;
360  }
361 }
std::map< dof_id_type, dof_id_type > _slave_to_master_map
maps slave node ids to master element ids
bool shouldApply() override
Whether or not this constraint should be applied.
SubProblem & _subproblem
SubProblem that contains tag info.
ConstraintType
Definition: MooseTypes.h:523
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:2267
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
Definition: FEProblem.h:24
DenseMatrix< Number > _Kee
stiffness matrix holding slave-slave jacobian
unsigned int number() const
Get variable number coming from libMesh.
void resize(unsigned int size)
Change the number of elements the array can store.
Definition: MooseArray.h:218
virtual Real computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, unsigned int jvar) override
This is the virtual that derived classes should override for computing the off-diag Jacobian...
InputParameters validParams< EqualValueEmbeddedConstraint >()
Real _constraint_residual
constraint force needed to enforce the constraint
T & set(const std::string &name, bool quiet_mode=false)
Returns a writable reference to the named parameters.
void reinitConstraint()
Prepare the residual contribution of the current constraint required to enforce it based on the speci...
const VariableValue & _u_master
Holds the current solution at the current quadrature point.
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
virtual void computeOffDiagJacobian(unsigned int jvar) override
Computes d-residual / d-jvar...
NumericVector< Number > & _residual_copy
copy of the residual before the constraint is applied
unsigned int _i
Definition: Constraint.h:75
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
Formulation
Formulations, currently only supports KINEMATIC and PENALTY.
virtual void getConnectedDofIndices(unsigned int var_num) override
Gets the indices for all dofs connected to the constraint.
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...
const VariableTestValue & _test_master
Side test function.
const std::string & type() const
Get the type of this object.
Definition: MooseObject.h:53
MooseVariableFEBase & getVariable(THREAD_ID tid, const std::string &var_name)
Gets a reference to a variable of with specified name.
Definition: SystemBase.C:105
virtual const dof_id_type & nodalDofIndex() const =0
MooseVariable & _var
virtual void getConnectedDofIndices(unsigned int var_num)
Gets the indices for all dofs connected to the constraint.
unsigned int _j
Definition: Constraint.h:75
unsigned short _master
master block id
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:259
static MooseEnum getNonlinearVariableOrders()
Get the possible variable orders.
registerMooseObject("MooseApp", EqualValueEmbeddedConstraint)
const Real _penalty
Penalty parameter used in constraint enforcement for kinematic and penalty formulations.
MeshBase & getMesh()
Accessor for the underlying libMesh Mesh object.
Definition: MooseMesh.C:2567
MooseMesh & _mesh
Definition: Constraint.h:73
const VariablePhiValue & _phi_master
Side shape function.
enum EqualValueEmbeddedConstraint::Formulation _formulation
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
VariablePhiValue _phi_slave
Shape function on the slave side.
unsigned short _slave
slave block id
virtual unsigned int number() const
Gets the number of this system.
Definition: SystemBase.C:926
EqualValueEmbeddedConstraint(const InputParameters &parameters)
virtual void computeJacobian() override
Computes the jacobian for the current element.
Assembly & _assembly
Definition: Constraint.h:72
const VariableValue & _u_slave
Value of the unknown variable on the slave node.
VariableTestValue _test_slave
Shape function on the slave side. This will always only have one entry and that entry will always be ...
virtual Real computeQpResidual(Moose::ConstraintType type) override
This is the virtual that derived classes should override for computing the residual.
MatType type
const Node *const & _current_node
current node being processed
virtual void prepareSlaveToMasterMap() override
prepare the _slave_to_master_map
virtual void addGhostedElem(dof_id_type elem_id)=0
Will make sure that all dofs connected to elem_id are ghosted to this processor.
SystemBase & _sys
Definition: Constraint.h:68
ConstraintJacobianType
Definition: MooseTypes.h:543
std::vector< dof_id_type > _connected_dof_indices
dofs connected to the slave node
virtual Real computeQpJacobian(Moose::ConstraintJacobianType type) override
This is the virtual that derived classes should override for computing the Jacobian.
MooseVariable & _master_var
Master side variable.
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...
DenseMatrix< Number > & jacobianBlockNeighbor(Moose::DGJacobianType type, unsigned int ivar, unsigned int jvar, TagID tag=0)
Definition: Assembly.C:2033
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
virtual Real computeQpSlaveValue() override
Compute the value the slave node should have at the beginning of a timestep.
A EqualValueEmbeddedConstraint forces the value of a variable to be the same on overlapping portion o...
A NodeElemConstraint is used when you need to create constraints between a slave node and a master el...
DenseMatrix< Number > _Kne
stiffness matrix holding master-slave jacobian
virtual std::unique_ptr< PointLocatorBase > getPointLocator() const
Proxy function to get a (sub)PointLocator from either the underlying libMesh mesh (default)...
Definition: MooseMesh.C:2730
InputParameters validParams< NodeElemConstraint >()
unsigned int _qp
Definition: Constraint.h:76
bool _overwrite_slave_residual
Whether or not the slave&#39;s residual should be overwritten.
virtual void setNeighborSubdomainID(const Elem *elem, unsigned int side, THREAD_ID tid) override
virtual void reinitNeighborPhys(const Elem *neighbor, unsigned int neighbor_side, const std::vector< Point > &physical_points, THREAD_ID tid) override