Adding a test case to show volumetric constraint using Lagrange multi…

…pliers (refs #000)

r23211

framework/src/base/Assembly.C

@@ -1172,6 +1172,7 @@ Assembly::addJacobian(SparseMatrix<Number> & jacobian)
         // We only add jacobian blocks for d(nl-var)/d(scalar-var) now (these we generated by kernels and BCs)
         // jacobian blocks d(scalar-var)/d(nl-var) are added later with addJacobianScalar
         addJacobianBlock(jacobian, jacobianBlock(jvar.index(), ivar.index()), jvar.dofIndices(), ivar.dofIndices(), jvar.scalingFactor());
+        addJacobianBlock(jacobian, jacobianBlock(ivar.index(), jvar.index()), ivar.dofIndices(), jvar.dofIndices(), ivar.scalingFactor());
       }
     }
   }

test/include/kernels/ScalarLagrangeMultiplier.h

@@ -0,0 +1,30 @@
+#ifndef SCALARLAGRANGEMULTIPLIER_H
+#define SCALARLAGRANGEMULTIPLIER_H
+
+#include "Kernel.h"
+
+class ScalarLagrangeMultiplier;
+
+template<>
+InputParameters validParams<ScalarLagrangeMultiplier>();
+
+/**
+ * Implements coupling of Lagrange multiplier (as a scalar variable) to a simple constraint of type (g(u) - c = 0)
+ */
+class ScalarLagrangeMultiplier : public Kernel
+{
+public:
+  ScalarLagrangeMultiplier(const std::string & name, InputParameters parameters);
+  virtual ~ScalarLagrangeMultiplier();
+
+  virtual void computeOffDiagJacobianScalar(unsigned int jvar);
+
+protected:
+  virtual Real computeQpResidual();
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar);
+
+  unsigned int _lambda_var;
+  VariableValue & _lambda;
+};
+
+#endif /* SCALARLAGRANGEMULTIPLIER_H */

test/include/scalarkernels/PostprocessorCED.h

@@ -0,0 +1,37 @@
+#ifndef POSTPROCESSORCED_H
+#define POSTPROCESSORCED_H
+
+#include "ScalarKernel.h"
+
+class PostprocessorCED;
+
+template<>
+InputParameters validParams<PostprocessorCED>();
+
+/**
+ *
+ */
+class PostprocessorCED : public ScalarKernel
+{
+public:
+  PostprocessorCED(const std::string & name, InputParameters parameters);
+  virtual ~PostprocessorCED();
+
+  virtual void reinit();
+  virtual void computeResidual();
+  virtual void computeJacobian();
+  virtual void computeOffDiagJacobian(unsigned int jvar);
+
+protected:
+  virtual Real computeQpResidual();
+  virtual Real computeQpJacobian();
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar);
+
+  unsigned int _i;
+
+  Real _value;
+  PostprocessorValue & _pp_value;
+};
+
+
+#endif /* POSTPROCESSORCED_H */

test/src/base/MooseTestApp.C

@@ -41,6 +41,7 @@
 #include "ConvectionPrecompute.h"
 #include "CoupledKernelGradTest.h"
 #include "CoupledKernelValueTest.h"
+#include "ScalarLagrangeMultiplier.h"
 #include "SplineFFn.h"
 #include "BlkResTestDiffusion.h"
 #include "DiffTensorKernel.h"
@@ -97,6 +98,7 @@
 #include "ImplicitODEx.h"
 #include "ImplicitODEy.h"
 #include "AlphaCED.h"
+#include "PostprocessorCED.h"
 
 #include "EqualValueNodalConstraint.h"
 // user objects
@@ -221,6 +223,7 @@ MooseTestApp::registerObjects(Factory & factory)
   registerKernel(SplineFFn);
   registerKernel(BlkResTestDiffusion);
   registerKernel(DiffTensorKernel);
+  registerKernel(ScalarLagrangeMultiplier);
 
   // Aux kernels
   registerAux(CoupledAux);
@@ -283,6 +286,7 @@ MooseTestApp::registerObjects(Factory & factory)
   registerScalarKernel(ImplicitODEx);
   registerScalarKernel(ImplicitODEy);
   registerScalarKernel(AlphaCED);
+  registerScalarKernel(PostprocessorCED);
 
   // Functions
   registerFunction(TimestepSetupFunction);

test/src/kernels/ScalarLagrangeMultiplier.C

@@ -0,0 +1,53 @@
+#include "ScalarLagrangeMultiplier.h"
+
+template<>
+InputParameters validParams<ScalarLagrangeMultiplier>()
+{
+  InputParameters params = validParams<Kernel>();
+  params.addRequiredCoupledVar("lambda", "Lagrange multiplier variable");
+
+  return params;
+}
+
+ScalarLagrangeMultiplier::ScalarLagrangeMultiplier(const std::string & name, InputParameters parameters) :
+    Kernel(name, parameters),
+    _lambda_var(coupledScalar("lambda")),
+    _lambda(coupledScalarValue("lambda"))
+{
+}
+
+ScalarLagrangeMultiplier::~ScalarLagrangeMultiplier()
+{
+}
+
+Real
+ScalarLagrangeMultiplier::computeQpResidual()
+{
+  return _lambda[0] * _test[_i][_qp];
+}
+
+void
+ScalarLagrangeMultiplier::computeOffDiagJacobianScalar(unsigned int jvar)
+{
+  DenseMatrix<Number> & ken = _assembly.jacobianBlock(_var.index(), jvar);
+  DenseMatrix<Number> & kne = _assembly.jacobianBlock(jvar, _var.index());
+  MooseVariableScalar & jv = _sys.getScalarVariable(_tid, jvar);
+
+  for (_i = 0; _i < _test.size(); _i++)
+    for (_j = 0; _j < jv.order(); _j++)
+      for (_qp = 0; _qp < _qrule->n_points(); _qp++)
+      {
+        Real value = _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar);
+        ken(_i, _j) += value;
+        kne(_j, _i) += value;
+      }
+}
+
+Real
+ScalarLagrangeMultiplier::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  if (jvar == _lambda_var)
+    return _test[_i][_qp];
+  else
+    return 0.;
+}

test/src/scalarkernels/PostprocessorCED.C

@@ -0,0 +1,66 @@
+#include "PostprocessorCED.h"
+
+template<>
+InputParameters validParams<PostprocessorCED>()
+{
+  InputParameters params = validParams<ScalarKernel>();
+  params.addRequiredParam<PostprocessorName>("pp_name", "");
+  params.addRequiredParam<Real>("value", "");
+
+  return params;
+}
+
+PostprocessorCED::PostprocessorCED(const std::string & name, InputParameters parameters) :
+    ScalarKernel(name, parameters),
+    _value(getParam<Real>("value")),
+    _pp_value(getPostprocessorValue("pp_name"))
+{
+}
+
+PostprocessorCED::~PostprocessorCED()
+{
+}
+
+void
+PostprocessorCED::reinit()
+{
+}
+
+void
+PostprocessorCED::computeResidual()
+{
+  DenseVector<Number> & re = _assembly.residualBlock(_var.index());
+  for (_i = 0; _i < re.size(); _i++)
+    re(_i) += computeQpResidual();
+}
+
+Real
+PostprocessorCED::computeQpResidual()
+{
+  return _pp_value - _value;
+}
+
+void
+PostprocessorCED::computeJacobian()
+{
+  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.index(), _var.index());
+  for (_i = 0; _i < ke.m(); _i++)
+    ke(_i, _i) += computeQpJacobian();
+}
+
+Real
+PostprocessorCED::computeQpJacobian()
+{
+  return 0.;
+}
+
+void
+PostprocessorCED::computeOffDiagJacobian(unsigned int /*jvar*/)
+{
+}
+
+Real
+PostprocessorCED::computeQpOffDiagJacobian(unsigned int /*jvar*/)
+{
+  return 0.;
+}

test/tests/kernels/scalar_constraint/scalar_constraint_kernel.i

@@ -0,0 +1,145 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  xmin = -1
+  xmax = 1
+  ymin = -1
+  ymax = 1
+  nx = 2
+  ny = 2
+  elem_type = QUAD9
+[]
+
+[Functions]
+  [./exact_fn]
+    type = ParsedFunction
+    value = 'x*x+y*y'
+  [../]
+  
+  [./ffn]
+    type = ParsedFunction
+    value = -4
+  [../]
+  
+  [./bottom_bc_fn]
+    type = ParsedFunction
+    value = -2*y
+  [../]
+
+  [./right_bc_fn]
+    type = ParsedFunction
+    value =  2*x
+  [../]
+
+  [./top_bc_fn]
+    type = ParsedFunction
+    value =  2*y
+  [../]
+
+  [./left_bc_fn]
+    type = ParsedFunction
+    value = -2*x
+  [../]
+[]
+
+# NL
+
+[Variables]
+  [./u]
+    family = LAGRANGE
+    order = SECOND
+  [../]
+
+  [./lambda]
+    family = SCALAR
+    order = FIRST
+  [../]
+[]
+
+[Kernels]
+  [./diff]
+    type = Diffusion
+    variable = u
+  [../]
+
+  [./ffnk]
+    type = UserForcingFunction
+    variable = u
+    function = ffn
+  [../]
+
+  [./sk_lm]
+    type = ScalarLagrangeMultiplier
+    variable = u
+    lambda = lambda
+  [../]
+[]
+
+[ScalarKernels]
+  [./constraint]
+    type = PostprocessorCED
+    variable = lambda
+    pp_name = pp
+    value = 2.666666666666666
+  [../]
+[]
+
+[BCs]
+  [./bottom]
+    type = FunctionNeumannBC
+    variable = u
+    boundary = '0'
+    function = bottom_bc_fn
+  [../]
+  [./right]
+    type = FunctionNeumannBC
+    variable = u
+    boundary = '1'
+    function = right_bc_fn
+  [../]
+  [./top]
+    type = FunctionNeumannBC
+    variable = u
+    boundary = '2'
+    function = top_bc_fn
+  [../]
+  [./left]
+    type = FunctionNeumannBC
+    variable = u
+    boundary = '3'
+    function = left_bc_fn
+  [../]
+[]
+
+[Postprocessors]
+  [./pp]
+    type = ElementIntegralVariablePostprocessor
+    variable = u
+    execute_on = residual
+  [../]
+  [./l2_err]
+    type = ElementL2Error
+    variable = u
+    function = exact_fn
+  [../]
+[]
+
+[Preconditioning]
+  active = 'pc'
+  
+  [./pc]
+    type = SMP
+    full = true
+    solve_type = 'PJFNK'
+  [../]
+[] # End preconditioning block
+
+[Executioner]
+  type = Steady
+  nl_rel_tol = 1e-15
+[]
+
+[Output]
+  output_initial = true
+  exodus = true
+[]

test/tests/kernels/scalar_constraint/tests

@@ -1,4 +1,11 @@
 [Tests]
+  [./kernel]
+    type = 'Exodiff'
+    input = 'scalar_constraint_kernel.i'
+    exodiff = 'scalar_constraint_kernel_out.e'
+    max_parallel = 1
+  [../]
+
   [./bc]
     type = 'Exodiff'
     input = 'scalar_constraint_bc.i'