add array variable action and array initial condition and a simple test

idaholab#6881

framework/include/actions/AddVariableAction.h

@@ -74,6 +74,9 @@ class AddVariableAction : public Action, public OutputInterface
   /// True if the variable being created is a scalar
   bool _scalar_var;
 
+  /// Number of components for an array variable
+  unsigned int _component;
+
   /// Absolute zero tolerance
   static const Real _abs_zero_tol;
 };

framework/include/ics/ArrayConstantIC.h

@@ -0,0 +1,36 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifndef ARRAYCONSTANTIC_H
+#define ARRAYCONSTANTIC_H
+
+#include "ArrayInitialCondition.h"
+
+// Forward Declarations
+class ArrayConstantIC;
+namespace libMesh
+{
+class Point;
+}
+
+template <>
+InputParameters validParams<ArrayConstantIC>();
+
+class ArrayConstantIC : public ArrayInitialCondition
+{
+public:
+  ArrayConstantIC(const InputParameters & parameters);
+
+  virtual RealArrayValue value(const Point & p) override;
+
+protected:
+  const RealArrayValue _value;
+};
+
+#endif // VECTORCONSTANTIC_H

framework/include/ics/ArrayInitialCondition.h

@@ -0,0 +1,15 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifndef ARRAYINITIALCONDITION_H
+#define ARRAYINITIALCONDITION_H
+
+#include "InitialConditionTempl.h"
+
+#endif // ARRAYINITIALCONDITION_H

framework/include/ics/InitialConditionTempl.h

@@ -26,6 +26,8 @@ template <>
 InputParameters validParams<InitialConditionTempl<Real>>();
 template <>
 InputParameters validParams<InitialConditionTempl<RealVectorValue>>();
+template <>
+InputParameters validParams<InitialConditionTempl<RealArrayValue>>();
 
 /**
  * This is a template class that implements the workhorse `compute` and `computeNodal` methods. The
@@ -37,9 +39,11 @@ template <typename T>
 class InitialConditionTempl : public InitialConditionBase
 {
 public:
-  typedef FEGenericBase<T> FEBaseType;
   typedef typename OutputTools<T>::OutputShape ValueType;
   typedef typename OutputTools<T>::OutputGradient GradientType;
+  typedef typename OutputTools<T>::OutputShapeGradient GradientShapeType;
+  typedef typename OutputTools<T>::OutputData DataType;
+  typedef FEGenericBase<ValueType> FEBaseType;
 
   /**
    * Constructor
@@ -60,7 +64,7 @@ class InitialConditionTempl : public InitialConditionBase
    *
    * This must be overridden by derived classes.
    */
-  virtual ValueType value(const Point & p) = 0;
+  virtual T value(const Point & p) = 0;
 
   /**
    * The gradient of the variable at a point.
@@ -70,6 +74,8 @@ class InitialConditionTempl : public InitialConditionBase
    */
   virtual GradientType gradient(const Point & /*p*/) { return GradientType(); };
 
+  T gradientComponent(GradientType grad, unsigned int i);
+
   /**
    * set the temporary solution vector for node projections of C0 variables
    */
@@ -87,13 +93,26 @@ class InitialConditionTempl : public InitialConditionBase
    * Helps perform multiplication of GradientTypes: a normal dot product for vectors and a
    * contraction for tensors
    */
-  Real dotHelper(const GradientType & op1, const GradientType & op2);
+  Real dotHelper(const RealGradient & op1, const RealGradient & op2) { return op1 * op2; }
+  Real dotHelper(const RealTensor & op1, const RealTensor & op2) { return op1.contract(op2); }
+  RealArrayValue dotHelper(const RealVectorArrayValue & op1, const RealGradient & op2)
+  {
+    RealArrayValue v = op1.col(0) * op2(0);
+    for (unsigned int i = 1; i < LIBMESH_DIM; ++i)
+      v += op1.col(i) * op2(i);
+    return v;
+  }
 
   /**
    * Perform the cholesky solves for edge, side, and interior projections
    */
   void choleskySolve(bool is_volume);
 
+  /**
+   * Assemble a small local system for cholesky solve
+   */
+  void choleskyAssembly(bool is_volume);
+
 protected:
   FEProblemBase & _fe_problem;
   THREAD_ID _tid;
@@ -123,9 +142,9 @@ class InitialConditionTempl : public InitialConditionBase
   /// Matrix storage member
   DenseMatrix<Real> _Ke;
   /// Linear b vector
-  DenseVector<Number> _Fe;
+  DenseVector<DataType> _Fe;
   /// Linear solution vector
-  DenseVector<Number> _Ue;
+  DenseVector<DataType> _Ue;
 
   /// The finite element type for the IC variable
   const FEType & _fe_type;
@@ -157,7 +176,7 @@ class InitialConditionTempl : public InitialConditionBase
   /// pointers to shape functions
   const std::vector<std::vector<ValueType>> * _phi;
   /// pointers to shape function gradients
-  const std::vector<std::vector<GradientType>> * _dphi;
+  const std::vector<std::vector<GradientShapeType>> * _dphi;
   /// pointers to the Jacobian * quadrature weights for current element
   const std::vector<Real> * _JxW;
   /// pointers to the xyz coordinates of the quadrature points for the current element
@@ -171,3 +190,4 @@ class InitialConditionTempl : public InitialConditionBase
 
 typedef InitialConditionTempl<Real> InitialCondition;
 typedef InitialConditionTempl<RealVectorValue> VectorInitialCondition;
+typedef InitialConditionTempl<RealArrayValue> ArrayInitialCondition;

framework/include/postprocessors/ElementIntegralArrayVariablePostprocessor.h

@@ -0,0 +1,45 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#ifndef ELEMENTINTEGRALVARIABLEARRAYPOSTPROCESSOR_H
+#define ELEMENTINTEGRALVARIABLEARRAYPOSTPROCESSOR_H
+
+#include "ElementIntegralPostprocessor.h"
+#include "MooseVariableInterface.h"
+
+// Forward Declarations
+class ElementIntegralArrayVariablePostprocessor;
+
+template <>
+InputParameters validParams<ElementIntegralArrayVariablePostprocessor>();
+
+/**
+ * This postprocessor computes a volume integral of the specified variable.
+ *
+ * Note that specializations of this integral are possible by deriving from this
+ * class and overriding computeQpIntegral().
+ */
+class ElementIntegralArrayVariablePostprocessor : public ElementIntegralPostprocessor,
+                                                  public MooseVariableInterface<RealArrayValue>
+{
+public:
+  ElementIntegralArrayVariablePostprocessor(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpIntegral() override;
+
+  /// Holds the solution at current quadrature points
+  const ArrayVariableValue & _u;
+  /// Holds the solution gradient at the current quadrature points
+  const ArrayVariableGradient & _grad_u;
+  /// The component
+  const unsigned int _component;
+};
+
+#endif

framework/src/actions/AddAuxVariableAction.C

@@ -29,7 +29,9 @@ validParams<AddAuxVariableAction>()
                              "Specifies the order of the FE shape function to use "
                              "for this variable (additional orders not listed are "
                              "allowed)");
-  params.addParam<Real>("initial_condition", "Specifies the initial condition for this variable");
+  params.addParam<unsigned int>("component", 1, "Number of component for an array variable");
+  params.addParam<std::vector<Real>>("initial_condition",
+                                     "Specifies the initial condition for this variable");
   params.addParam<std::vector<SubdomainName>>("block", "The block id where this variable lives");
 
   return params;

framework/src/actions/AddVariableAction.C

@@ -50,12 +50,15 @@ validParams<AddVariableAction>()
                              "Specifies the order of the FE shape function to use "
                              "for this variable (additional orders not listed are "
                              "allowed)");
-  params.addParam<Real>("initial_condition", "Specifies the initial condition for this variable");
+  params.addParam<unsigned int>("component", 1, "Number of component for an array variable");
+  params.addParam<std::vector<Real>>("initial_condition",
+                                     "Specifies the initial condition for this variable");
   params.addParam<std::vector<SubdomainName>>("block", "The block id where this variable lives");
   params.addParam<bool>("eigen", false, "True to make this variable an eigen variable");
 
   // Advanced input options
-  params.addParam<Real>("scaling", 1.0, "Specifies a scaling factor to apply to this variable");
+  params.addParam<std::vector<Real>>("scaling",
+                                     "Specifies a scaling factor to apply to this variable");
   params.addParamNamesToGroup("scaling eigen", "Advanced");
 
   return params;
@@ -66,7 +69,8 @@ AddVariableAction::AddVariableAction(InputParameters params)
     OutputInterface(params, false),
     _fe_type(Utility::string_to_enum<Order>(getParam<MooseEnum>("order")),
              Utility::string_to_enum<FEFamily>(getParam<MooseEnum>("family"))),
-    _scalar_var(_fe_type.family == SCALAR)
+    _scalar_var(_fe_type.family == SCALAR),
+    _component(getParam<unsigned int>("component"))
 {
 }
 
@@ -113,16 +117,29 @@ AddVariableAction::createInitialConditionAction()
 
   if (_scalar_var)
     action_params.set<std::string>("type") = "ScalarConstantIC";
-  else
+  else if (_component == 1)
     action_params.set<std::string>("type") = "ConstantIC";
+  else
+    action_params.set<std::string>("type") = "ArrayConstantIC";
 
   // Create the action
   std::shared_ptr<MooseObjectAction> action = std::static_pointer_cast<MooseObjectAction>(
       _action_factory.create("AddInitialConditionAction", long_name, action_params));
 
   // Set the required parameters for the object to be created
   action->getObjectParams().set<VariableName>("variable") = var_name;
-  action->getObjectParams().set<Real>("value") = getParam<Real>("initial_condition");
+  auto value = getParam<std::vector<Real>>("initial_condition");
+  if (value.size() != _component)
+    mooseError("Size of 'initial_condition' is not consistent");
+  if (_component > 1)
+  {
+    RealArrayValue v(_component);
+    for (unsigned int i = 0; i < _component; ++i)
+      v(i) = value[i];
+    action->getObjectParams().set<RealArrayValue>("value") = v;
+  }
+  else
+    action->getObjectParams().set<Real>("value") = value[0];
 
   // Store the action in the ActionWarehouse
   _awh.addActionBlock(action);
@@ -132,19 +149,34 @@ void
 AddVariableAction::addVariable(const std::string & var_name)
 {
   std::set<SubdomainID> blocks = getSubdomainIDs();
-  Real scale_factor = isParamValid("scaling") ? getParam<Real>("scaling") : 1;
+  std::vector<Real> scale_factor = isParamValid("scaling") ? getParam<std::vector<Real>>("scaling")
+                                                           : std::vector<Real>(_component, 1);
+  if (scale_factor.size() != _component)
+    mooseError("Size of 'scaling' is not consistent");
 
   // Scalar variable
   if (_scalar_var)
-    _problem->addScalarVariable(var_name, _fe_type.order, scale_factor);
+    _problem->addScalarVariable(var_name, _fe_type.order, scale_factor[0]);
 
   // Block restricted variable
-  else if (blocks.empty())
-    _problem->addVariable(var_name, _fe_type, scale_factor);
+  else if (_component == 1)
+  {
+    if (blocks.empty())
+      _problem->addVariable(var_name, _fe_type, scale_factor[0]);
 
-  // Non-block restricted variable
+    // Non-block restricted variable
+    else
+      _problem->addVariable(var_name, _fe_type, scale_factor[0], &blocks);
+  }
   else
-    _problem->addVariable(var_name, _fe_type, scale_factor, &blocks);
+  {
+    if (blocks.empty())
+      _problem->addArrayVariable(var_name, _fe_type, _component, scale_factor);
+
+    // Non-block restricted variable
+    else
+      _problem->addArrayVariable(var_name, _fe_type, _component, scale_factor, &blocks);
+  }
 
   if (getParam<bool>("eigen"))
   {

framework/src/ics/ArrayConstantIC.C

@@ -0,0 +1,37 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ArrayConstantIC.h"
+
+#include "libmesh/point.h"
+
+registerMooseObject("MooseApp", ArrayConstantIC);
+
+template <>
+InputParameters
+validParams<ArrayConstantIC>()
+{
+  InputParameters params = validParams<ArrayInitialCondition>();
+  params.addRequiredParam<RealArrayValue>("value", "The values to be set in IC");
+  params.addClassDescription("Sets constant component values for an array field variable.");
+  return params;
+}
+
+ArrayConstantIC::ArrayConstantIC(const InputParameters & parameters)
+  : ArrayInitialCondition(parameters), _value(getParam<RealArrayValue>("value"))
+{
+  if (_var.count() != _value.size())
+    mooseError("'value' size is inconsistent to the number of components of the array variable");
+}
+
+RealArrayValue
+ArrayConstantIC::value(const Point & /*p*/)
+{
+  return _value;
+}