Line data    Source code

       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             : #include "NodalConstraint.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "MooseVariableFE.h"
      15             : #include "SystemBase.h"
      16             : 
      17             : #include "libmesh/sparse_matrix.h"
      18             : 
      19             : InputParameters
      20        9346 : NodalConstraint::validParams()
      21             : {
      22        9346 :   InputParameters params = Constraint::validParams();
      23       37384 :   MooseEnum formulationtype("penalty kinematic", "penalty");
      24       37384 :   params.addParam<MooseEnum>("formulation",
      25             :                              formulationtype,
      26             :                              "Formulation used to calculate constraint - penalty or kinematic.");
      27       28038 :   params.addParam<NonlinearVariableName>("variable_secondary",
      28             :                                          "The name of the variable for the secondary nodes, if it "
      29             :                                          "is different from the primary nodes' variable");
      30       18692 :   return params;
      31        9346 : }
      32             : 
      33          88 : NodalConstraint::NodalConstraint(const InputParameters & parameters)
      34             :   : Constraint(parameters),
      35             :     NeighborCoupleableMooseVariableDependencyIntermediateInterface(this, true, true),
      36             :     NeighborMooseVariableInterface<Real>(
      37             :         this, true, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_STANDARD),
      38          88 :     _var(_sys.getFieldVariable<Real>(_tid, parameters.get<NonlinearVariableName>("variable"))),
      39          88 :     _var_secondary(_sys.getFieldVariable<Real>(
      40             :         _tid,
      41         176 :         isParamValid("variable_secondary")
      42          88 :             ? parameters.get<NonlinearVariableName>("variable_secondary")
      43         169 :             : parameters.get<NonlinearVariableName>("variable"))),
      44          88 :     _u_secondary(_var_secondary.dofValuesNeighbor()),
      45         176 :     _u_primary(_var.dofValues())
      46             : {
      47          88 :   addMooseVariableDependency(&_var);
      48          88 :   addMooseVariableDependency(&_var_secondary);
      49             : 
      50         176 :   MooseEnum temp_formulation = getParam<MooseEnum>("formulation");
      51          88 :   if (temp_formulation == "penalty")
      52          88 :     _formulation = Moose::Penalty;
      53           0 :   else if (temp_formulation == "kinematic")
      54           0 :     _formulation = Moose::Kinematic;
      55             :   else
      56           0 :     mooseError("Formulation must be either Penalty or Kinematic");
      57          88 : }
      58             : 
      59             : void
      60         850 : NodalConstraint::computeResidual(const NumericVector<Number> & residual)
      61             : {
      62         850 :   if ((_weights.size() == 0) && (_primary_node_vector.size() == 1))
      63          50 :     _weights.push_back(1.0);
      64             : 
      65         850 :   std::vector<dof_id_type> primarydof = _var.dofIndices();
      66         850 :   std::vector<dof_id_type> secondarydof = _var_secondary.dofIndicesNeighbor();
      67             : 
      68         850 :   DenseVector<Number> re(primarydof.size());
      69         850 :   DenseVector<Number> neighbor_re(secondarydof.size());
      70             : 
      71         850 :   re.zero();
      72         850 :   neighbor_re.zero();
      73             : 
      74        3858 :   for (_i = 0; _i < secondarydof.size(); ++_i)
      75             :   {
      76        6016 :     for (_j = 0; _j < primarydof.size(); ++_j)
      77             :     {
      78        3008 :       switch (_formulation)
      79             :       {
      80        3008 :         case Moose::Penalty:
      81        3008 :           re(_j) += computeQpResidual(Moose::Primary) * _var.scalingFactor();
      82        3008 :           neighbor_re(_i) += computeQpResidual(Moose::Secondary) * _var_secondary.scalingFactor();
      83        3008 :           break;
      84           0 :         case Moose::Kinematic:
      85             :           // Transfer the current residual of the secondary node to the primary nodes
      86           0 :           Real res = residual(secondarydof[_i]);
      87           0 :           re(_j) += res * _weights[_j];
      88           0 :           neighbor_re(_i) +=
      89           0 :               -res / _primary_node_vector.size() + computeQpResidual(Moose::Secondary);
      90           0 :           break;
      91             :       }
      92             :     }
      93             :   }
      94             :   // We've already applied scaling
      95         850 :   if (!primarydof.empty())
      96         826 :     addResiduals(_assembly, re, primarydof, /*scaling_factor=*/1);
      97         850 :   if (!secondarydof.empty())
      98         826 :     addResiduals(_assembly, neighbor_re, secondarydof, /*scaling_factor=*/1);
      99         850 : }
     100             : 
     101             : void
     102         149 : NodalConstraint::computeJacobian(const SparseMatrix<Number> & jacobian)
     103             : {
     104         149 :   if ((_weights.size() == 0) && (_primary_node_vector.size() == 1))
     105           0 :     _weights.push_back(1.0);
     106             : 
     107             :   // Calculate the dense-block Jacobian entries
     108         149 :   std::vector<dof_id_type> secondarydof = _var_secondary.dofIndicesNeighbor();
     109         149 :   std::vector<dof_id_type> primarydof = _var.dofIndices();
     110             : 
     111         149 :   DenseMatrix<Number> Kee(primarydof.size(), primarydof.size());
     112         149 :   DenseMatrix<Number> Ken(primarydof.size(), secondarydof.size());
     113         149 :   DenseMatrix<Number> Kne(secondarydof.size(), primarydof.size());
     114             : 
     115         149 :   Kee.zero();
     116         149 :   Ken.zero();
     117         149 :   Kne.zero();
     118             : 
     119         792 :   for (_i = 0; _i < secondarydof.size(); ++_i)
     120             :   {
     121        1286 :     for (_j = 0; _j < primarydof.size(); ++_j)
     122             :     {
     123         643 :       switch (_formulation)
     124             :       {
     125         643 :         case Moose::Penalty:
     126         643 :           Kee(_j, _j) += computeQpJacobian(Moose::PrimaryPrimary);
     127         643 :           Ken(_j, _i) += computeQpJacobian(Moose::PrimarySecondary);
     128         643 :           Kne(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
     129         643 :           break;
     130           0 :         case Moose::Kinematic:
     131           0 :           Kee(_j, _j) = 0.;
     132           0 :           Ken(_j, _i) += jacobian(secondarydof[_i], primarydof[_j]) * _weights[_j];
     133           0 :           Kne(_i, _j) += -jacobian(secondarydof[_i], primarydof[_j]) / primarydof.size() +
     134           0 :                          computeQpJacobian(Moose::SecondaryPrimary);
     135           0 :           break;
     136             :       }
     137             :     }
     138             :   }
     139         149 :   addJacobian(_assembly, Kee, primarydof, primarydof, _var.scalingFactor());
     140         149 :   addJacobian(_assembly, Ken, primarydof, secondarydof, _var.scalingFactor());
     141         149 :   addJacobian(_assembly, Kne, secondarydof, primarydof, _var_secondary.scalingFactor());
     142             : 
     143             :   // Calculate and cache the diagonal secondary-secondary entries
     144         792 :   for (_i = 0; _i < secondarydof.size(); ++_i)
     145             :   {
     146         643 :     Number value = 0.0;
     147         643 :     switch (_formulation)
     148             :     {
     149         643 :       case Moose::Penalty:
     150         643 :         value = computeQpJacobian(Moose::SecondarySecondary);
     151         643 :         break;
     152           0 :       case Moose::Kinematic:
     153           0 :         value = -jacobian(secondarydof[_i], secondarydof[_i]) / primarydof.size() +
     154           0 :                 computeQpJacobian(Moose::SecondarySecondary);
     155           0 :         break;
     156             :     }
     157         643 :     addJacobianElement(
     158         643 :         _assembly, value, secondarydof[_i], secondarydof[_i], _var_secondary.scalingFactor());
     159             :   }
     160         149 : }
     161             : 
     162             : void
     163           0 : NodalConstraint::updateConnectivity()
     164             : {
     165           0 : }