From e9f751b77f2c9f66102883089c202ffab9d8246b Mon Sep 17 00:00:00 2001 From: Yaqi Wang Date: Wed, 24 Apr 2019 11:03:11 -0600 Subject: [PATCH] add array variable action and array initial condition and a simple test #6881 --- framework/include/actions/AddVariableAction.h | 3 + framework/include/ics/ArrayConstantIC.h | 36 ++++ framework/include/ics/ArrayInitialCondition.h | 15 ++ framework/include/ics/InitialConditionTempl.h | 32 +++- ...lementIntegralArrayVariablePostprocessor.h | 45 +++++ framework/src/actions/AddAuxVariableAction.C | 4 +- framework/src/actions/AddVariableAction.C | 54 ++++-- framework/src/ics/ArrayConstantIC.C | 37 ++++ framework/src/ics/InitialConditionTempl.C | 174 +++++++++++------- ...lementIntegralArrayVariablePostprocessor.C | 40 ++++ framework/src/problems/FEProblemBase.C | 2 + .../array_variable/array_variable_test.i | 141 ++++++++++++++ .../gold/array_variable_test_out.e | Bin 0 -> 77740 bytes test/tests/variables/array_variable/tests | 11 ++ 14 files changed, 512 insertions(+), 82 deletions(-) create mode 100644 framework/include/ics/ArrayConstantIC.h create mode 100644 framework/include/ics/ArrayInitialCondition.h create mode 100644 framework/include/postprocessors/ElementIntegralArrayVariablePostprocessor.h create mode 100644 framework/src/ics/ArrayConstantIC.C create mode 100644 framework/src/postprocessors/ElementIntegralArrayVariablePostprocessor.C create mode 100644 test/tests/variables/array_variable/array_variable_test.i create mode 100644 test/tests/variables/array_variable/gold/array_variable_test_out.e create mode 100644 test/tests/variables/array_variable/tests diff --git a/framework/include/actions/AddVariableAction.h b/framework/include/actions/AddVariableAction.h index cdafa82968a6..d93c3a485eba 100644 --- a/framework/include/actions/AddVariableAction.h +++ b/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; }; diff --git a/framework/include/ics/ArrayConstantIC.h b/framework/include/ics/ArrayConstantIC.h new file mode 100644 index 000000000000..7273a5355f19 --- /dev/null +++ b/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(); + +class ArrayConstantIC : public ArrayInitialCondition +{ +public: + ArrayConstantIC(const InputParameters & parameters); + + virtual RealArrayValue value(const Point & p) override; + +protected: + const RealArrayValue _value; +}; + +#endif // VECTORCONSTANTIC_H diff --git a/framework/include/ics/ArrayInitialCondition.h b/framework/include/ics/ArrayInitialCondition.h new file mode 100644 index 000000000000..f2d387ac6288 --- /dev/null +++ b/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 diff --git a/framework/include/ics/InitialConditionTempl.h b/framework/include/ics/InitialConditionTempl.h index 0c06e84107ea..ce5187e0d3e8 100644 --- a/framework/include/ics/InitialConditionTempl.h +++ b/framework/include/ics/InitialConditionTempl.h @@ -26,6 +26,8 @@ template <> InputParameters validParams>(); template <> InputParameters validParams>(); +template <> +InputParameters validParams>(); /** * This is a template class that implements the workhorse `compute` and `computeNodal` methods. The @@ -37,9 +39,11 @@ template class InitialConditionTempl : public InitialConditionBase { public: - typedef FEGenericBase FEBaseType; typedef typename OutputTools::OutputShape ValueType; typedef typename OutputTools::OutputGradient GradientType; + typedef typename OutputTools::OutputShapeGradient GradientShapeType; + typedef typename OutputTools::OutputData DataType; + typedef FEGenericBase 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 _Ke; /// Linear b vector - DenseVector _Fe; + DenseVector _Fe; /// Linear solution vector - DenseVector _Ue; + DenseVector _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> * _phi; /// pointers to shape function gradients - const std::vector> * _dphi; + const std::vector> * _dphi; /// pointers to the Jacobian * quadrature weights for current element const std::vector * _JxW; /// pointers to the xyz coordinates of the quadrature points for the current element @@ -171,3 +190,4 @@ class InitialConditionTempl : public InitialConditionBase typedef InitialConditionTempl InitialCondition; typedef InitialConditionTempl VectorInitialCondition; +typedef InitialConditionTempl ArrayInitialCondition; diff --git a/framework/include/postprocessors/ElementIntegralArrayVariablePostprocessor.h b/framework/include/postprocessors/ElementIntegralArrayVariablePostprocessor.h new file mode 100644 index 000000000000..a8e8c2660006 --- /dev/null +++ b/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(); + +/** + * 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 +{ +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 diff --git a/framework/src/actions/AddAuxVariableAction.C b/framework/src/actions/AddAuxVariableAction.C index f0292e7cc9a1..3a421ef4d599 100644 --- a/framework/src/actions/AddAuxVariableAction.C +++ b/framework/src/actions/AddAuxVariableAction.C @@ -29,7 +29,9 @@ validParams() "Specifies the order of the FE shape function to use " "for this variable (additional orders not listed are " "allowed)"); - params.addParam("initial_condition", "Specifies the initial condition for this variable"); + params.addParam("component", 1, "Number of component for an array variable"); + params.addParam>("initial_condition", + "Specifies the initial condition for this variable"); params.addParam>("block", "The block id where this variable lives"); return params; diff --git a/framework/src/actions/AddVariableAction.C b/framework/src/actions/AddVariableAction.C index 903d1ed6a188..2c59c8916f40 100644 --- a/framework/src/actions/AddVariableAction.C +++ b/framework/src/actions/AddVariableAction.C @@ -50,12 +50,15 @@ validParams() "Specifies the order of the FE shape function to use " "for this variable (additional orders not listed are " "allowed)"); - params.addParam("initial_condition", "Specifies the initial condition for this variable"); + params.addParam("component", 1, "Number of component for an array variable"); + params.addParam>("initial_condition", + "Specifies the initial condition for this variable"); params.addParam>("block", "The block id where this variable lives"); params.addParam("eigen", false, "True to make this variable an eigen variable"); // Advanced input options - params.addParam("scaling", 1.0, "Specifies a scaling factor to apply to this variable"); + params.addParam>("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(getParam("order")), Utility::string_to_enum(getParam("family"))), - _scalar_var(_fe_type.family == SCALAR) + _scalar_var(_fe_type.family == SCALAR), + _component(getParam("component")) { } @@ -113,8 +117,10 @@ AddVariableAction::createInitialConditionAction() if (_scalar_var) action_params.set("type") = "ScalarConstantIC"; - else + else if (_component == 1) action_params.set("type") = "ConstantIC"; + else + action_params.set("type") = "ArrayConstantIC"; // Create the action std::shared_ptr action = std::static_pointer_cast( @@ -122,7 +128,18 @@ AddVariableAction::createInitialConditionAction() // Set the required parameters for the object to be created action->getObjectParams().set("variable") = var_name; - action->getObjectParams().set("value") = getParam("initial_condition"); + auto value = getParam>("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("value") = v; + } + else + action->getObjectParams().set("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 blocks = getSubdomainIDs(); - Real scale_factor = isParamValid("scaling") ? getParam("scaling") : 1; + std::vector scale_factor = isParamValid("scaling") ? getParam>("scaling") + : std::vector(_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("eigen")) { diff --git a/framework/src/ics/ArrayConstantIC.C b/framework/src/ics/ArrayConstantIC.C new file mode 100644 index 000000000000..e934fc886b62 --- /dev/null +++ b/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() +{ + InputParameters params = validParams(); + params.addRequiredParam("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("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; +} diff --git a/framework/src/ics/InitialConditionTempl.C b/framework/src/ics/InitialConditionTempl.C index 26a96015438e..49b116d127e2 100644 --- a/framework/src/ics/InitialConditionTempl.C +++ b/framework/src/ics/InitialConditionTempl.C @@ -28,6 +28,13 @@ validParams>() return validParams(); } +template <> +InputParameters +validParams>() +{ + return validParams(); +} + template InitialConditionTempl::InitialConditionTempl(const InputParameters & parameters) : InitialConditionBase(parameters), @@ -90,7 +97,7 @@ InitialConditionTempl::compute() if (_cont == C_ONE) { - const std::vector> & ref_dphi = fe->get_dphi(); + const std::vector> & ref_dphi = fe->get_dphi(); _dphi = &ref_dphi; } @@ -228,31 +235,12 @@ InitialConditionTempl::compute() NumericVector & solution = _var.sys().solution(); - // 'first' and 'last' are no longer used, see note about subdomain-restricted variables below - // const dof_id_type - // first = solution.first_local_index(), - // last = solution.last_local_index(); - // Lock the new_vector since it is shared among threads. { Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx); - for (unsigned int i = 0; i < n_dofs; i++) - // We may be projecting a new zero value onto - // an old nonzero approximation - RHS - // if (_Ue(i) != 0.) - - // This is commented out because of subdomain restricted variables. - // It can be the case that if a subdomain restricted variable's boundary - // aligns perfectly with a processor boundary that the variable will get - // no value. To counteract this we're going to let every processor set a - // value at every node and then let PETSc figure it out. - // Later we can choose to do something different / better. - // if ((_dof_indices[i] >= first) && (_dof_indices[i] < last)) - { - solution.set(_dof_indices[i], _Ue(i)); - } _var.setDofValues(_Ue); + _var.insert(solution); } } @@ -285,6 +273,27 @@ InitialConditionTempl::setCZeroVertices() } } +template +T +InitialConditionTempl::gradientComponent(GradientType grad, unsigned int i) +{ + return grad(i); +} + +template <> +RealVectorValue +InitialConditionTempl::gradientComponent(GradientType grad, unsigned int i) +{ + return grad.row(i); +} + +template <> +RealArrayValue +InitialConditionTempl::gradientComponent(GradientType grad, unsigned int i) +{ + return grad.col(i); +} + template void InitialConditionTempl::setHermiteVertices() @@ -295,9 +304,9 @@ InitialConditionTempl::setHermiteVertices() _Ue(_current_dof) = value(*_current_node); _dof_is_fixed[_current_dof] = true; _current_dof++; - Gradient grad = gradient(*_current_node); + GradientType grad = gradient(*_current_node); // x derivative - _Ue(_current_dof) = grad(0); + _Ue(_current_dof) = gradientComponent(grad, 0); _dof_is_fixed[_current_dof] = true; _current_dof++; if (_dim > 1) @@ -306,35 +315,38 @@ InitialConditionTempl::setHermiteVertices() Point nxminus = _current_elem->point(_n), nxplus = _current_elem->point(_n); nxminus(0) -= TOLERANCE; nxplus(0) += TOLERANCE; - Gradient gxminus = gradient(nxminus); - Gradient gxplus = gradient(nxplus); + GradientType gxminus = gradient(nxminus); + GradientType gxplus = gradient(nxplus); // y derivative - _Ue(_current_dof) = grad(1); + _Ue(_current_dof) = gradientComponent(grad, 1); _dof_is_fixed[_current_dof] = true; _current_dof++; // xy derivative - _Ue(_current_dof) = (gxplus(1) - gxminus(1)) / 2. / TOLERANCE; + _Ue(_current_dof) = + (gradientComponent(gxplus, 1) - gradientComponent(gxminus, 1)) / 2. / TOLERANCE; _dof_is_fixed[_current_dof] = true; _current_dof++; if (_dim > 2) { // z derivative - _Ue(_current_dof) = grad(2); + _Ue(_current_dof) = gradientComponent(grad, 2); _dof_is_fixed[_current_dof] = true; _current_dof++; // xz derivative - _Ue(_current_dof) = (gxplus(2) - gxminus(2)) / 2. / TOLERANCE; + _Ue(_current_dof) = + (gradientComponent(gxplus, 2) - gradientComponent(gxminus, 2)) / 2. / TOLERANCE; _dof_is_fixed[_current_dof] = true; _current_dof++; // We need new points for yz Point nyminus = _current_elem->point(_n), nyplus = _current_elem->point(_n); nyminus(1) -= TOLERANCE; nyplus(1) += TOLERANCE; - Gradient gyminus = gradient(nyminus); - Gradient gyplus = gradient(nyplus); + GradientType gyminus = gradient(nyminus); + GradientType gyplus = gradient(nyplus); // xz derivative - _Ue(_current_dof) = (gyplus(2) - gyminus(2)) / 2. / TOLERANCE; + _Ue(_current_dof) = + (gradientComponent(gyplus, 2) - gradientComponent(gyminus, 2)) / 2. / TOLERANCE; _dof_is_fixed[_current_dof] = true; _current_dof++; // Getting a 2nd order xyz is more tedious @@ -348,12 +360,14 @@ InitialConditionTempl::setHermiteVertices() nxpym(1) -= TOLERANCE; nxpyp(0) += TOLERANCE; nxpyp(1) += TOLERANCE; - Gradient gxmym = gradient(nxmym); - Gradient gxmyp = gradient(nxmyp); - Gradient gxpym = gradient(nxpym); - Gradient gxpyp = gradient(nxpyp); - Number gxzplus = (gxpyp(2) - gxmyp(2)) / 2. / TOLERANCE; - Number gxzminus = (gxpym(2) - gxmym(2)) / 2. / TOLERANCE; + GradientType gxmym = gradient(nxmym); + GradientType gxmyp = gradient(nxmyp); + GradientType gxpym = gradient(nxpym); + GradientType gxpyp = gradient(nxpyp); + DataType gxzplus = + (gradientComponent(gxpyp, 2) - gradientComponent(gxmyp, 2)) / 2. / TOLERANCE; + DataType gxzminus = + (gradientComponent(gxpym, 2) - gradientComponent(gxmym, 2)) / 2. / TOLERANCE; // xyz derivative _Ue(_current_dof) = (gxzplus - gxzminus) / 2. / TOLERANCE; _dof_is_fixed[_current_dof] = true; @@ -380,10 +394,10 @@ InitialConditionTempl::setOtherCOneVertices() _Ue(_current_dof) = value(*_current_node); _dof_is_fixed[_current_dof] = true; _current_dof++; - Gradient grad = gradient(*_current_node); + GradientType grad = gradient(*_current_node); for (unsigned int i = 0; i != _dim; ++i) { - _Ue(_current_dof) = grad(i); + _Ue(_current_dof) = gradientComponent(grad, i); _dof_is_fixed[_current_dof] = true; _current_dof++; } @@ -395,32 +409,10 @@ InitialConditionTempl::setOtherCOneVertices() { } -template -Real -InitialConditionTempl::dotHelper(const GradientType & op1, const GradientType & op2) -{ - return op1 * op2; -} - -template <> -Real -InitialConditionTempl::dotHelper(const GradientType & op1, - const GradientType & op2) -{ - return op1.contract(op2); -} - template void -InitialConditionTempl::choleskySolve(bool is_volume) +InitialConditionTempl::choleskyAssembly(bool is_volume) { - _Ke.resize(_free_dofs, _free_dofs); - _Ke.zero(); - _Fe.resize(_free_dofs); - _Fe.zero(); - // The new edge coefficients - DenseVector U(_free_dofs); - // Loop over the quadrature points for (_qp = 0; _qp < _n_qp; _qp++) { @@ -464,6 +456,21 @@ InitialConditionTempl::choleskySolve(bool is_volume) freei++; } } +} + +template +void +InitialConditionTempl::choleskySolve(bool is_volume) +{ + _Ke.resize(_free_dofs, _free_dofs); + _Ke.zero(); + _Fe.resize(_free_dofs); + _Fe.zero(); + + choleskyAssembly(is_volume); + + // The new edge coefficients + DenseVector U(_free_dofs); _Ke.cholesky_solve(_Fe, U); @@ -471,13 +478,51 @@ InitialConditionTempl::choleskySolve(bool is_volume) for (unsigned int i = 0; i != _free_dofs; ++i) { auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]]; - Number & ui = _Ue(the_dof); + DataType & ui = _Ue(the_dof); libmesh_assert(std::abs(ui) < TOLERANCE || std::abs(ui - U(i)) < TOLERANCE); ui = U(i); _dof_is_fixed[the_dof] = true; } } +template <> +void +InitialConditionTempl::choleskySolve(bool is_volume) +{ + _Ke.resize(_free_dofs, _free_dofs); + _Ke.zero(); + _Fe.resize(_free_dofs); + for (unsigned int i = 0; i < _free_dofs; ++i) + _Fe(i).setZero(_var.count()); + + choleskyAssembly(is_volume); + + // The new edge coefficients + DenseVector U = _Fe; + + for (unsigned int i = 0; i < _var.count(); ++i) + { + DenseVector v(_free_dofs), x(_free_dofs); + for (unsigned int j = 0; j < _free_dofs; ++j) + v(j) = _Fe(j)(i); + + _Ke.cholesky_solve(v, x); + + for (unsigned int j = 0; j < _free_dofs; ++j) + U(j)(i) = x(j); + } + + // Transfer new edge solutions to element + for (unsigned int i = 0; i != _free_dofs; ++i) + { + auto the_dof = is_volume ? _free_dof[i] : _side_dofs[_free_dof[i]]; + DataType & ui = _Ue(the_dof); + libmesh_assert(ui.matrix().norm() < TOLERANCE || (ui - U(i)).matrix().norm() < TOLERANCE); + ui = U(i); + _dof_is_fixed[the_dof] = true; + } +} + template void InitialConditionTempl::computeNodal(const Point & p) @@ -498,3 +543,4 @@ InitialConditionTempl::computeNodal(const Point & p) template class InitialConditionTempl; template class InitialConditionTempl; +template class InitialConditionTempl; diff --git a/framework/src/postprocessors/ElementIntegralArrayVariablePostprocessor.C b/framework/src/postprocessors/ElementIntegralArrayVariablePostprocessor.C new file mode 100644 index 000000000000..3f094e33f7c4 --- /dev/null +++ b/framework/src/postprocessors/ElementIntegralArrayVariablePostprocessor.C @@ -0,0 +1,40 @@ +//* 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 "ElementIntegralArrayVariablePostprocessor.h" + +registerMooseObject("MooseApp", ElementIntegralArrayVariablePostprocessor); + +template <> +InputParameters +validParams() +{ + InputParameters params = validParams(); + params.addRequiredCoupledVar("variable", "The name of the variable that this object operates on"); + params.addParam("component", 0, ""); + return params; +} + +ElementIntegralArrayVariablePostprocessor::ElementIntegralArrayVariablePostprocessor( + const InputParameters & parameters) + : ElementIntegralPostprocessor(parameters), + MooseVariableInterface( + this, false, "variable", Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_ARRAY), + _u(coupledArrayValue("variable")), + _grad_u(coupledArrayGradient("variable")), + _component(getParam("component")) +{ + addMooseVariableDependency(mooseVariable()); +} + +Real +ElementIntegralArrayVariablePostprocessor::computeQpIntegral() +{ + return _u[_qp](_component); +} diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index 068e66e22138..8a7808d76c26 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -2256,6 +2256,8 @@ FEProblemBase::addInitialCondition(const std::string & ic_name, ic = _factory.create(ic_name, name, parameters, tid); else if (dynamic_cast(&var)) ic = _factory.create(ic_name, name, parameters, tid); + else if (dynamic_cast(&var)) + ic = _factory.create(ic_name, name, parameters, tid); else mooseError("Your FE variable in initial condition ", name, diff --git a/test/tests/variables/array_variable/array_variable_test.i b/test/tests/variables/array_variable/array_variable_test.i new file mode 100644 index 000000000000..c0c77693fe95 --- /dev/null +++ b/test/tests/variables/array_variable/array_variable_test.i @@ -0,0 +1,141 @@ +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 10 + ny = 8 +[] + +[Problem] + kernel_coverage_check = false + solve = false +[] + +[Variables] + [./u] + order = FIRST + family = LAGRANGE + component = 4 + initial_condition = '1 2 3 4' + [../] + [./uu] + order = FIRST + family = LAGRANGE + component = 2 + initial_condition = '1 2' + [../] + [./v] + order = FIRST + family = LAGRANGE + component = 2 + initial_condition = '5 6' + [../] + [./w] + order = CONSTANT + family = MONOMIAL + component = 3 + initial_condition = '7 8 9' + [../] + [./x] + order = THIRD + family = MONOMIAL + component = 2 + initial_condition = '10 11' + [../] + [./y] + order = FIRST + family = L2_LAGRANGE + component = 3 + initial_condition = '12 13 14' + [../] +[] + +[Postprocessors] + [u0int] + type = ElementIntegralArrayVariablePostprocessor + variable = u + component = 0 + [] + [u1int] + type = ElementIntegralArrayVariablePostprocessor + variable = u + component = 1 + [] + [u2int] + type = ElementIntegralArrayVariablePostprocessor + variable = u + component = 2 + [] + [u3int] + type = ElementIntegralArrayVariablePostprocessor + variable = u + component = 3 + [] + [uu0int] + type = ElementIntegralArrayVariablePostprocessor + variable = uu + component = 0 + [] + [uu1int] + type = ElementIntegralArrayVariablePostprocessor + variable = uu + component = 1 + [] + [v0int] + type = ElementIntegralArrayVariablePostprocessor + variable = v + component = 0 + [] + [v1int] + type = ElementIntegralArrayVariablePostprocessor + variable = v + component = 1 + [] + [w0int] + type = ElementIntegralArrayVariablePostprocessor + variable = w + component = 0 + [] + [w1int] + type = ElementIntegralArrayVariablePostprocessor + variable = w + component = 1 + [] + [w2int] + type = ElementIntegralArrayVariablePostprocessor + variable = w + component = 2 + [] + [x0int] + type = ElementIntegralArrayVariablePostprocessor + variable = x + component = 0 + [] + [x1int] + type = ElementIntegralArrayVariablePostprocessor + variable = x + component = 1 + [] + [y0int] + type = ElementIntegralArrayVariablePostprocessor + variable = y + component = 0 + [] + [y1int] + type = ElementIntegralArrayVariablePostprocessor + variable = y + component = 1 + [] + [y2int] + type = ElementIntegralArrayVariablePostprocessor + variable = y + component = 2 + [] +[] + +[Executioner] + type = Steady +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/variables/array_variable/gold/array_variable_test_out.e b/test/tests/variables/array_variable/gold/array_variable_test_out.e new file mode 100644 index 0000000000000000000000000000000000000000..9f579ef1ded428922a088c6e18fce31609309fbd GIT binary patch literal 77740 zcmeHQOK=>=d0qf~fHWyl7H!McYbeSR^`h_+CCY{|Fds5(J|q#6Ez|aDu{!|QTeP&K)aS(y*4>$U$?#*phB>QgbAvAv$D$atiQL}TP@e>tIX!H_A%XT)q)NsrBW3% zPFc_hSm;&Iew{8aHLfC6Dvf{+O1*L?=$g9JLCHJy8uW*vfWUTuMIZji-l|B)d5w5f^E+d>RBV4@hwF=E& zItWtw(Qaie9AjA?ccsxXqS7JAlzr>XlNS1bkM*Kh7 zKMTG$MgG8YqgCm`L&|Hddb3-;-fGuiw|Cgzv?Ke9cG|aGZ=g?fMqeKD51`_S88}X* z+hbQsJw!{pyY+5^U9l(~PkvR+@MPr>KPaMqRX z71k9@4{_R_2H&sz_4<{TIRJH`%_QtWgt{aA^jZAd7c50eA}wuW0R|p}gmCaM{*9yj zTD!GI|B|4WJD${6 z?duF4kG?CR|9D@7a2)0PkEg~6Lt||XRmPfn$%d_1n^i6QF%JI-`~c^~G6(N!WesV&&8Ry<-4PDo#J}-stJw^yUHTf^OdmgB zh;{d=d7-| zAQarZvg<5saA{WJ8h=Mj_qVab!+(IZ;O3S63Cn7K2<>|xWD=_iy3b2kL2u$*K3NMbtaMiC_F0Y1a7^i2g4DvXLL^#3SDaKd9=ciNUPa zbiGd_&#G6d>3SZ!>ON&lsW%xcS>uRvyabG?(eLTR|24qlb*GBf{mm~2>P{D}dy~i2 z=M6No+^sAjP5X?4*j5Q_Gy5FD`(ymmHjNgktn(R=Wt&kmmwJE7I?6WjI$XQ4o&?S_ z-_wcPY!}v<%dK)dsAAR8L7K;T0K|D(0%Skv>cnF|0+&4Ueqs?Gk22!)=P}OEpRWXS zKmJ38KN7Bl_&4_#MlSqsr&X8}{`V$d_|J9CIpP0e$_xM3WfkUx|A(nB{Fi^p@KX$3o84xc?czMA~g_l5F!oLjq0*LEh-qZUMXb$9q=0RTu zodTT(odKN%y#o3QXaV$9(APlcK<7ahKo>!;f?fl?4!Q(d1ib-z6ZCb^H$ZQJz6rVv z`WEQhpznanptnI4&=RN$8Uc-h#z5=~#-%?fV_f=0`bhdm`bzG9NS{f6N#9AoNgqo8 zNnc7oN}o!9vTvyq^^m@0oO-a2DPx>^u)ir|oO-bDDPx>^&<`kMoO;kFC}W&@&_5_+ zoO;k#C}W&@cp%Ccr>`9Wc_8{5^>_?))WoR={faWisRwc4M?L3Gw>O|cx zfG(Oib$ivw7^iNp85!f$?R6t#oVr~yGRCRfqLDF9J>D=f#;M1fM#ebx__~oXPTe>U zQAVAp+gqS-nmBd4Y-Eg6w{IC4Tpd+BiKu19`7N0=;NzgISQ=q3o&w!o##pitF&4a!SIt8LHiNAW4~}5ar}sX9YLJ_z= zmO#~Bc>*jW^vd+FjL<97z!{-eX76Q$UY-&tkgdeEyjEPx0Vu8KA*|;iY~&$qYZu_m2Ob;mNvZe7cX8q;~mX7J^?fFN(l;bCe5vv ziVHyRRBLs$(yV#suocIfZ?E)LgJ!qmA?prsz2nE9*MF<6R!48Dc&@e9^`5JH3HvRM z*QvpqdZUQ9fSc;}anH8#D+#uOTvfU4RqPOHHDAuLQ)+urPxpg4PzQ?&J}4z1X-z8r>0(-IdO7ga(bpZ z{ZeK6r58_3OjT5o& zKGkct(Wl-bHdRIHq2-HLdfxn6+naj9o1A%h=H$!MCsb?o!kqPe{zLH9f8l4=gK7`I z$}|FniS0J4YYNx`oOpS9S_8)G3(hv;^wWjrTCeM!#r7lPr5`E%FLNQZm3H>Ldb3jP z)~^N8Ac@6hy<^Q*xfU!}dW|lZvN^{mV7*(fG|J2E)@oVp$`&16WSO5DsuYwSK%{`C5kz%Zr@xxxw2PBZQ-uxSj zSgo4!x!m(ofGohGc)sVUEidV)TEoJ4juf8)}nn4PlPt^F`=*R1mW z_WxkYUDop#P9#B_g9f3vy^_igKT>jz2f)t7rSD=F);tCv&OvsV_4Iil2PbS-fnL7U z4Jx$_>9<3sooU?jaC!^YFi@Jj==hf;RwVBm)M`*%k ztG!yTUJ0sKwe7NTLj5rxHq~i9K{+J~{cxZ^ll{l@8s*9oH-`pm<;Zdo{H*sJJ-eWO zC7*!wzdHIcF+jr!go66{hX|OYSO^Zl;JFE`1sMvg!vKhtwei0Yr7tq;nNw&nyiyFX zyjH1Rt*pc*@L6vSJJzcO^(z687F=J(cJJO=4HqnZ%BkqAH#2o|=H!VNr%o1XHvpHH zYvuG4%dB^Ojp`M#AB%ZjGuv-er@WAJ|Df-9Az)fp8^8knOZ`f73Xo=8?lmz%4F=dG@R~M0?Q=1=tQ6sbf!UsM=FN&k!)vh7-I`JGk!BW{Ar@?YCuIFXEh+ zV8m#G^^@uk9LEC0nFSY<{aq6R7Q=%L)yNVJu_Lw7O9m{YA06Aw)SMKP>}WcrpnjYW z4IoaLRc-vjR{_nSv$~8eh&|k8Ndh$M!{qt{xBzT5uHk-7+OF(USZ}V>vE(M;5O^hC z_zY8?ntZV3FY`u~y5iZQ1#lb0h~}zYFKUb@JLOgRmH^#9s^vt(&|KxuKh49pWMzqeB2`F zb#dWi!9+?bKP-+{^1LOecC}ouMHsPI*eQed@6#*l9|5$_o*`%im{@p^_Dgc{QB_xg`3^WS6mK`VX8e>@2g&bZvrm*0Iu4@mc!`lzRH?ereVi`;XN`ocv| zTR){YP65~;->jwGOH=mK`ixupCAS-ZSLDzo*y%8g3%I=N}~t`3;l4Y02Ni zZ5H0J0305mTMGC9bsSPpU+h1Iy%rxeY0+=5kC9TZbjs^-l)jV#l=6Zr{ViPD4f`6v zjXsEYlz`HoC8LknPx?bRk81$miKFzT066Y*7YX(gSp0xV?>Bsm9Ag=Gbla_ZEpi<^ zqFt-w`3-e*EgE|Fn0mrTO5R+6Ef-D{^ymS-qJV5)0c0?>EUtOMr)~p{mkW+wsrpFX>^aNziH}=e$wuK zyP64BJ&<1M;?^M_Iqx@QL$^=+0$P+1;;(y(xu3GLTXPk1seF2J$w>eS;W zl0FsDcn-&LP`oy3CUs3>}UK;1NNsm}~ z@WW>M#9%Z8pwusg^~DYv|5kWBT#Z)Rm}7Y@nt1Iq>zzJxcK(fXi*Y?|drJF{DL{RG zO8NnbHjHlfNq>kv?e?-A1Wa;2wpIn)otp^TBYN0e>pQ}70%$IOu*$*tZKbcDX5m+O zNdp_4ire3Q&s)dM0i#V4gWDwQBkdPHs{4z8yaG!N<}|u_)95Gmlm5Q#Kc%|seV@Av z+2oB}4O!u)sVDU}OH8?`o--m`f~-_uy2zJ7`a-#_WMUNj_7+^FNyK`n|_e~xCWHt zOTD!I=k+pwy~UC^ZE=GpbDb&uE?n9j?{^#zfSqgo8rg9Ij`-Zg&UR}q=_HOdc54~0 zRIzcwGytVu5xf$X#A;tvfZx52yaKEbz>bDZ*5#A>A=d68?{IpAMgOLA0%F$%eu}bHGieUJjQRukMJ_9b z!7GJ#IP1&mr;4)Lm{E^@;R|h;)y_UWv4RcS%ecg=ti^BUiUmwBXWFd~4;rrF){5RK zf)`H6wuxxunE9}vejF(pkhhz$@1oLN!JR+8G$O8HVLV>|Rl^%C>n;#bRg;Io^Ey{r z*F~KK9(D>1iR|S2)NO#>Yrs2^kb9l z{fAllu~WK1HD_`dLO%jt>t0P%%?QO03$+*a#f49;QUoxiAEf=^@d#-B8XAW;8WzR# z_4-ioOa2kHp4_zROTOr1b30!-*B1S|{3?Ep!OP^TL8sGdt7YOAA#W9o=UqPjyn67% z>^}T$f!U>*dvym#O7& zJ=UcPm38#3@fyze*nEhMgL?VI2?&2?7onW+BUMsj1XRD9F_IKw^I-}|znB=@(uiHA z*n6K4!pWh(Rz_R{X51&9^-yOUndT*J=B7Or{vaXy4&ud;vfR7O&` z6AksxtUe5f<`c434A4{T=vtBR1C;Cco`a1sn7EN7!xLX(1 z6D%JH+afJ#H#fm}7M;A$do(n_zQINgKKemR5+u$900`1GR2~m44ol zD_D8nK6hSF1!yWa!HUCZ?pCn5zC%Op$Q3N#9^ocf8816>g3Woqj5T&>DmTH3!)Wd% zSj$NUlQh(hoM5-`CM<2}ZV_z1TLfx9P!RCx+dWYMl6@RnH zU0A_>EGJm?_Hj4C?&t~D-Gg-#tZe}fh!1)n>~&j-%=9){!K$w}xCwSgPq1zUn>mu* zg9dnm-@v=^?!pT8<2k{)eXugQlJCe{1a2R!n_z7Va4@(PtZkLy(&|Yr_4P2ff~7w( z?8ph$-6F^tIP?s65B9F#gI&+*hIsZixkaG9zUU@c4k3meIl;PH1o5FekfLX}3D!-p z1CEg`P%BvV^;$QZ-2}@a#IPeL zShs?W58Z(jJ;P0~Zh{?fjBJ5k!MY#0ID91R$O+c{$dx@%O7spl!G7E)Sn~?p9AEe_ zvoF8eZY|+;s#~-m=JuG(CvUmnv|q)X-F^UOyHVm~ zzZ|uVPpRTt@9nZi?x?S%Nte)u3m4v;KeupNjFiK3UeS9efSGS|AmFf$U2V168Q(si z^}g}eH*=cVpN-E<1k06PBlY{B!*A2deCIElJnYa5+U*v)xLN*grP{*x1)KQpbn~5{ z-NKvn`I2TlSBbYW1Gv~dFn+fJAEn|uILr8SO1q9%c}G5j*zj%>0bYZH}KIK(e3GfvHd_pGD zba+5?1wKHDPu?a&4(~S20N+|@_8Ryk{TenP;){%{L9;_oQDJxu)3NRlF#I;$Ql;Il zV;aNv96A+z9a?{XVR#L5zXLw(eAs>m-N6Z`vYzxEz*+Ataj`<&d}6gytEn%nqO(O~ zd>n~iCK$fS>YJ@i^zjBx1go|i-U%STRIXhl>HqK=Vy`;ZbfNkGq=B$@_VKMybi-4HA)LE*jCo$2E99IF8 z-jp}(&E&2BwdSmY*){n~o*MA%!t0k7v)kF1s{r)}>=+#k!>?^HXX8Wka!wjRyc0mN zq59yt`sP;`=Vr{`{Ms3Nc@bY9>uXhJlo|~J;&oV}>hKzFj}H)R;dDb&ebBSb@VEQi zY$Iz|?b>>_DK^t^2(gx%Z5o08=9imoGN;ujQ>|^TZS%E_n{A>FO^k7Jkkcq`S|zS+ zUhq!1*(T9?TG($xx30 zxixY(uB2Y{Uh+;Boo}dN&f4bG#S51f=Pwkzuuy>RN-8LS{I8_exA|=I(O=tagV|F@$!Tx0*lpqh zYa4f$xNo(Z6xF?~Z{0_iW*?=s&E5PNZth)e6UUt$XaVO3+7EoqN{%~w&;rg6j310U zP7C zlfT<=slP9!zTk+Swtk|w;Er3zm-=!%J74Oj;qEecf-iW&({M#!X`k?M|L7*eZNoPW zzTgNK`-)vfF8bK%lEJ%a`1^*xXYm_WxhP-zgY<)Nv6mcA>>z$1eAKi<^b@(1zgzY_ zV2*3!i+_v#q(6jKeBcbT;C117!S@G-+j4c~301viaLJ>h#yec^V! zgC<|(qL=W!*6~HZeMT<4WYTuNv|sRqOM6F5J&}tZ!tHj5T=bT9iJf-)pI77Jmf^Q8 zy~d1u)5vYT<#?jEw133Zd*ASzhTpLGZ>(~X|DHLX;0oVs%8Nc?4{6sPD_``Je6gF& zx5~a^AK{|6q(vX;7vUloJw<-de^T|&O~W?__%StZMt#joDLl%z@sjJOm)G?pJ(BB- z9;2bYqMxpR+vqtG=4<&aBe#!#+bSpZOQBr!h~wCLL^yGM?Dj~#h)cVro^X*17kQjZ zd+hc`?a<|;{CNFAJ>+=8rM_^H3l}^|3(u2FzHqyqEl3JGh*~kON)LHm*WY3JS})) zmy*F3dk*FcP9B%zXfFLBxN;ovlSrP7yW8j$`Lp2a@?xhb-^S79Z}`uc^ZB!8TxdSx z$J3)GEpp){BOf(5n}*-cslO}Cm*Z=`8Oq~&jD^RE%cJ9p+{Tf#97mV4+av8Lh3ybL z&GYhi6>PWQMD`Q?;>StL7rE#!cw&dh?izP69Y3DnNxN@_?T^=w%86aFac+k=y1eK+ zh$HQe<4L=s^6`9;3$AeS6R9ul5<5sak>9fXD$0LIjfW9m^OE6`&f`(}2i5qxZPkzR zHi}t<`@4ws6p>6N%;OsyP*nZ$^ zvy>dyhH!RpcAy2^esSXmkDTGWdY|8+7aBod1nV&!1;mY2X4Q-GKhO_!1#sp2kj5;`2pI&JwKqIXL!%-!w2<3-90}Lx1)Q0K=y&lK5t3u{oQx^{6e%J zJdP*(#iRZ5`Z(_S0ohkC_OiL;+gz6$NqH_~)5xReEu_BSOIoT+2cHj!+jB7A-j|=2mg5SZa4B!+i#(aT z=Lck8y1lPFo|b*$@?3+xUtIQ~+xyMsIfkP7?)icEa{%u70eLRK?f2LBF~lCl`1afH zA8%Kre_SrlGuS*TFaF{wv&p{nX#aj(F8MZ(%Desk+W(5#P3-9GC;Qyp^8?y{v^`|M zyUhpozlYCRIQ!j^{pS8#J-71fIsN#`ug$*9z5cyMZs*^x3!2|g_L}F*_W6%0`Hk>7 Z7=0ggGkgz)zRuihp2z)P^E~1B{{yL; literal 0 HcmV?d00001 diff --git a/test/tests/variables/array_variable/tests b/test/tests/variables/array_variable/tests new file mode 100644 index 000000000000..a66bb1fc5b78 --- /dev/null +++ b/test/tests/variables/array_variable/tests @@ -0,0 +1,11 @@ +[Tests] + [./block_aux_kernel_test] + type = 'Exodiff' + input = 'array_variable_test.i' + exodiff = 'array_variable_test_out.e' + group = 'requirements' + requirement = 'MOOSE shall provide an ability to add array variables with constant initial conditions.' + issues = '#6881' + design = 'moose_variables.md' + [../] +[]