diff --git a/framework/include/bcs/PenaltyDirichletBC.h b/framework/include/bcs/PenaltyDirichletBC.h new file mode 100644 index 000000000000..9b154c7bc604 --- /dev/null +++ b/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(); + +/** + * 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 & 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 _local_ke; +}; + +#endif /* PENALTYDIRICHLETBC_H */ diff --git a/framework/src/base/Moose.C b/framework/src/base/Moose.C index 8155229e17f4..a2cf50f283be 100644 --- a/framework/src/base/Moose.C +++ b/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); diff --git a/framework/src/bcs/PenaltyDirichletBC.C b/framework/src/bcs/PenaltyDirichletBC.C new file mode 100644 index 000000000000..988e46c68d01 --- /dev/null +++ b/framework/src/bcs/PenaltyDirichletBC.C @@ -0,0 +1,78 @@ +/****************************************************************/ +/* 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() +{ + InputParameters p = validParams(); + p.addRequiredParam("value", "Value of the BC"); + p.addRequiredParam("penalty", "Penalty factor"); + return p; +} + + +PenaltyDirichletBC::PenaltyDirichletBC(const std::string & name, InputParameters parameters) : + NodalBC(name, parameters), + _value(getParam("value")), + _penalty(getParam("penalty")) +{} + +void +PenaltyDirichletBC::computeResidual(NumericVector & 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() +{ + if (_var.isNodalDefined()) + { + _qp = 0; + Real jacobian_value = computeQpJacobian(); + dof_id_type dof_index = _var.nodalDofIndex(); + + DenseMatrix & ke = _assembly.jacobianBlock(_var.number(), _var.number()); + _local_ke.resize(ke.m(), ke.n()); + _local_ke.zero(); + + _local_ke(0,0) += computeQpJacobian(); + + ke += _local_ke; + } +} + +Real +PenaltyDirichletBC::computeQpJacobian() +{ + return _penalty; +} + +void +PenaltyDirichletBC::computeOffDiagJacobian() +{ +} diff --git a/modules/solid_mechanics/tests/penalty_dirichlet/gold/out.e b/modules/solid_mechanics/tests/penalty_dirichlet/gold/out.e new file mode 100644 index 000000000000..4ca9b47de030 Binary files /dev/null and b/modules/solid_mechanics/tests/penalty_dirichlet/gold/out.e differ diff --git a/modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i b/modules/solid_mechanics/tests/penalty_dirichlet/penalty_dirichlet.i new file mode 100644 index 000000000000..36781081591f --- /dev/null +++ b/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 + [../] +[] diff --git a/modules/solid_mechanics/tests/penalty_dirichlet/tests b/modules/solid_mechanics/tests/penalty_dirichlet/tests new file mode 100644 index 000000000000..f215c10afc67 --- /dev/null +++ b/modules/solid_mechanics/tests/penalty_dirichlet/tests @@ -0,0 +1,7 @@ +[Tests] + [./test] + type = 'Exodiff' + input = 'penalty_dirichlet.i' + exodiff = 'out.e' + [../] +[]