Add PenaltyDirichletBC closes idaholab#5268

framework/include/bcs/PenaltyDirichletBC.h

@@ -0,0 +1,52 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#ifndef PENALTYDIRICHLETBC_H
+#define PENALTYDIRICHLETBC_H
+
+#include "NodalBC.h"
+
+class PenaltyDirichletBC;
+
+template<>
+InputParameters validParams<PenaltyDirichletBC>();
+
+/**
+ * Boundary condition of a Dirichlet type
+ *
+ * Enforce using a penalty method
+ */
+class PenaltyDirichletBC : public NodalBC
+{
+public:
+  PenaltyDirichletBC(const std::string & name, InputParameters parameters);
+
+protected:
+  virtual void computeResidual(NumericVector<Number> & residual);
+  virtual Real computeQpResidual();
+  virtual void computeJacobian();
+  virtual Real computeQpJacobian();
+  virtual void computeOffDiagJacobian();
+
+  /// The value for this BC
+  const Real & _value;
+
+  /// The penalty factor
+  const Real & _penalty;
+
+  // Holds Jacobian entries
+  DenseMatrix<Number> _local_ke;
+};
+
+#endif /* PENALTYDIRICHLETBC_H */

framework/src/base/Moose.C

@@ -58,6 +58,7 @@
 // bcs
 #include "ConvectiveFluxBC.h"
 #include "DirichletBC.h"
+#include "PenaltyDirichletBC.h"
 #include "PresetBC.h"
 #include "NeumannBC.h"
 #include "FunctionDirichletBC.h"
@@ -427,6 +428,7 @@ registerObjects(Factory & factory)
   // bcs
   registerBoundaryCondition(ConvectiveFluxBC);
   registerBoundaryCondition(DirichletBC);
+  registerBoundaryCondition(PenaltyDirichletBC);
   registerBoundaryCondition(PresetBC);
   registerBoundaryCondition(NeumannBC);
   registerBoundaryCondition(FunctionDirichletBC);

framework/src/bcs/PenaltyDirichletBC.C

@@ -0,0 +1,85 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "PenaltyDirichletBC.h"
+
+template<>
+InputParameters validParams<PenaltyDirichletBC>()
+{
+  InputParameters p = validParams<NodalBC>();
+  p.addRequiredParam<Real>("value", "Value of the BC");
+  p.addRequiredParam<Real>("penalty", "Penalty factor");
+  return p;
+}
+
+
+PenaltyDirichletBC::PenaltyDirichletBC(const std::string & name, InputParameters parameters) :
+  NodalBC(name, parameters),
+  _value(getParam<Real>("value")),
+  _penalty(getParam<Real>("penalty"))
+{}
+
+void
+PenaltyDirichletBC::computeResidual(NumericVector<Number> & residual)
+{
+  if (_var.isNodalDefined())
+  {
+    dof_id_type & dof_idx = _var.nodalDofIndex();
+    _qp = 0;
+    residual.add(dof_idx, computeQpResidual());
+  }
+}
+
+Real
+PenaltyDirichletBC::computeQpResidual()
+{
+  return _penalty * (_u[_qp] - _value);
+}
+
+void
+PenaltyDirichletBC::computeJacobian()
+{
+  // We call the user's computeQpJacobian() function and store the
+  // results in the _assembly object. We can't store them directly in
+  // the element stiffness matrix, as they will only be inserted after
+  // all the assembly is done.
+  if (_var.isNodalDefined())
+  {
+    _qp = 0;
+    Real jacobian_value = computeQpJacobian();
+    dof_id_type dof_index = _var.nodalDofIndex();
+
+    DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), _var.number());
+    _local_ke.resize(ke.m(), ke.n());
+    _local_ke.zero();
+
+//    for (_qp = 0; _qp < _qrule->n_points(); _qp++)
+//      for (_i = 0; _i < _test.size(); _i++)
+//        for (_j = 0; _j < _phi.size(); _j++)
+//          _local_ke(_i, _j) += _JxW[_qp]*_coord[_qp]*computeQpJacobian();
+    _local_ke(0,0) += computeQpJacobian();
+
+    ke += _local_ke;
+  }
+}
+
+Real
+PenaltyDirichletBC::computeQpJacobian()
+{
+  return _penalty;
+}
+
+void
+PenaltyDirichletBC::computeOffDiagJacobian()
+{}

modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i

@@ -0,0 +1,124 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = 0.0
+  xmax = 1.0
+  ymin = 0.0
+  ymax = 1.0
+  nx = 1
+  ny = 1
+  displacements = 'disp_x disp_y'
+[]
+
+[Variables]
+  [./disp_x]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+  [./disp_y]
+    order = FIRST
+    family = LAGRANGE
+  [../]
+[]
+
+[AuxVariables]
+  [./stress_xx]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_yy]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./stress_zz]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+[]
+
+[SolidMechanics]
+  [./solid]
+    disp_x = disp_x
+    disp_y = disp_y
+  [../]
+[]
+
+[AuxKernels]
+  [./stress_xx]
+    type = MaterialTensorAux
+    variable = stress_xx
+    tensor = stress
+    index = 0
+  [../]
+  [./stress_yy]
+    type = MaterialTensorAux
+    variable = stress_yy
+    tensor = stress
+    index = 1
+  [../]
+  [./stress_zz]
+    type = MaterialTensorAux
+    variable = stress_zz
+    tensor = stress
+    index = 2
+  [../]
+[]
+
+[BCs]
+  [./left_x]
+    type = PenaltyDirichletBC
+    variable = disp_x
+    boundary = left
+    value = 0.0
+    penalty = 1e4
+  [../]
+  [./right_x]
+    type = PresetBC
+    variable = disp_x
+    boundary = right
+    value = 0.001
+  [../]
+  [./bottom_y]
+    type = DirichletBC
+    variable = disp_y
+    boundary = bottom
+    value = 0.0
+  [../]
+[]
+
+[Materials]
+  [./solid]
+    type = Elastic
+    block = 0
+    disp_x = disp_x
+    disp_y = disp_y
+    youngs_modulus = 2e4
+    poissons_ratio = 0.0
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+
+  #Preconditioned JFNK (default)
+  solve_type = 'PJFNK'
+
+  petsc_options = '-snes_ksp_ew '
+  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
+  petsc_options_value = '101 hypre boomeramg'
+  nl_rel_tol = 1e-12
+  l_tol = 1e-5
+  start_time = 0.0
+  dt = 1
+  num_steps = 1
+[]
+
+[Outputs]
+  file_base = out
+  output_initial = true
+  print_linear_residuals = true
+  print_perf_log = true
+  [./exodus]
+    type = Exodus
+  [../]
+[]

modules/solid_mechanics/tests/penalty_dirichlet/tests

@@ -0,0 +1,7 @@
+[Tests]
+  [./test]
+    type = 'Exodiff'
+    input = 'penalty_dirichlet.i'
+    exodiff = 'out.e'
+  [../]
+[]