Daniel Schwen's work on giving user objects access to shape functions (…

…idaholab#5913)

I'm trying to make Daniel's changes work on top of the warehouse changes.

framework/include/base/ComputeUserObjectsThread.h

@@ -60,6 +60,8 @@ class ComputeUserObjectsThread : public ThreadedElementLoop<ConstElemRange>
   const MooseObjectWarehouse<SideUserObject> & _side_user_objects;
   const MooseObjectWarehouse<InternalSideUserObject> & _internal_side_user_objects;
   ///@}
+
+  std::vector<MooseVariable *> _jacobian_moose_vars;
 };
 
 #endif //COMPUTEUSEROBJECTSTHREAD_H

framework/include/base/FEProblem.h

@@ -340,7 +340,7 @@ class FEProblem :
   /**
    * Returns true if we are currently computing Jacobian
    */
-  virtual bool currentlyComputingJacobian() { return _currently_computing_jacobian; }
+  virtual const bool & currentlyComputingJacobian() { return _currently_computing_jacobian; }
 
   /**
    * Returns true if we are in or beyond the initialSetup stage

framework/include/userobject/ShapeElementUserObject.h

@@ -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               */
+/****************************************************************/
+
+#ifndef SHAPEELEMENTUSEROBJECT_H
+#define SHAPEELEMENTUSEROBJECT_H
+
+#include "ElementUserObject.h"
+
+//Forward Declarations
+class ShapeElementUserObject;
+
+template<>
+InputParameters validParams<ShapeElementUserObject>();
+
+/**
+ * ElementUserObject class in which the _phi and _grad_phi shape function data
+ * is available and correctly initialized on EXEC_NONLINEAR (the Jacobian calculation).
+ * This enables the calculation of Jacobian matrix contributions inside a UO.
+ *
+ * \warning It is up to the user to ensure _fe_problem.currentlyComputingJacobian()
+ *          returns true before utilizing the shape functions.
+ */
+class ShapeElementUserObject : public ElementUserObject
+{
+public:
+  ShapeElementUserObject(const InputParameters & parameters);
+
+  /**
+   * Returns the set of variables a Jacobian has been requested for
+   */
+  const std::set<MooseVariable *> & jacobianMooseVariables() { return _jacobian_moose_variables; }
+
+  /**
+   * This function will be called with the shape functions for jvar initialized. It
+   * can be used to compute Jacobian contributions of the by implementing executeJacobian.
+   */
+  virtual void executeJacobianWrapper(unsigned int jvar, const std::vector<dof_id_type> & dof_indices);
+
+protected:
+  /**
+   * Implement this function to compute Jacobian terms for this UserObject. The
+   * shape function index _j and its corrsponding global DOF index _j_global
+   * will be provided.
+   */
+  virtual void executeJacobian(unsigned int /*jvar*/) = 0;
+
+  /**
+   * Returns the index for a coupled variable by name and requests the computation
+   * of a Jacobian w.r.t. to this variable i.e. the call to executeJacobian() with
+   * shapefunctions initialized for this variable.
+   */
+  virtual unsigned int coupled(const std::string & var_name, unsigned int comp = 0);
+
+  /// shape function values
+  const VariablePhiValue & _phi;
+
+  /// shape function gradients
+  const VariablePhiGradient & _grad_phi;
+
+  /// j-th index for enumerating the shape functions
+  unsigned int _j;
+
+  /// global DOF ID corresponding to _j
+  dof_id_type _j_global;
+
+  /// set to true iff the current call of the user object is for the purpose of calculating Jacobians
+  const bool & _currently_computing_jacobian;
+
+private:
+  const bool _compute_jacobians;
+  std::set<MooseVariable *> _jacobian_moose_variables;
+};
+
+#endif //SHAPEELEMENTUSEROBJECT_H

framework/src/base/ComputeUserObjectsThread.C

@@ -16,6 +16,7 @@
 #include "Problem.h"
 #include "SystemBase.h"
 #include "ElementUserObject.h"
+#include "ShapeElementUserObject.h"
 #include "SideUserObject.h"
 #include "InternalSideUserObject.h"
 #include "NodalUserObject.h"
@@ -64,6 +65,24 @@ ComputeUserObjectsThread::subdomainChanged()
 
   _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, _tid);
   _fe_problem.prepareMaterials(_subdomain, _tid);
+
+  // ShapeElementUserObject requested Jacobian variables
+  {
+    std::set<MooseVariable *> jacobian_moose_vars;
+
+    const std::vector<MooseSharedPointer<ElementUserObject> > & e_objects = _elemental_user_objects.getActiveBlockObjects(_subdomain, _tid);
+    for (std::vector<MooseSharedPointer<ElementUserObject> >::const_iterator it = e_objects.begin(); it != e_objects.end(); ++it)
+    {
+      MooseSharedPointer<ShapeElementUserObject> shape_element_uo = MooseSharedNamespace::dynamic_pointer_cast<ShapeElementUserObject>(*it);
+      if (shape_element_uo)
+      {
+        const std::set<MooseVariable *> & mv_deps = shape_element_uo->jacobianMooseVariables();
+        jacobian_moose_vars.insert(mv_deps.begin(), mv_deps.end());
+      }
+    }
+
+    _jacobian_moose_vars.assign(jacobian_moose_vars.begin(), jacobian_moose_vars.end());
+  }
 }
 
 void
@@ -80,6 +99,29 @@ ComputeUserObjectsThread::onElement(const Elem * elem)
       uo->execute();
   }
 
+  // UserObject Jacobians
+  if (_fe_problem.currentlyComputingJacobian())
+  {
+    // Prepare shape functions for ShapeElementUserObjects
+    for (std::vector<MooseVariable *>::const_iterator jvar_it = _jacobian_moose_vars.begin();
+         jvar_it != _jacobian_moose_vars.end();
+         ++jvar_it)
+    {
+      unsigned int jvar = (*jvar_it)->number();
+      std::vector<dof_id_type> & dof_indices = (*jvar_it)->dofIndices();
+
+      _fe_problem.prepareShapes(jvar, _tid);
+
+      const std::vector<MooseSharedPointer<ElementUserObject> > & e_objects = _elemental_user_objects.getActiveBlockObjects(_subdomain, _tid);
+      for (std::vector<MooseSharedPointer<ElementUserObject> >::const_iterator it = e_objects.begin(); it != e_objects.end(); ++it)
+      {
+        MooseSharedPointer<ShapeElementUserObject> shape_element_uo = MooseSharedNamespace::dynamic_pointer_cast<ShapeElementUserObject>(*it);
+        if (shape_element_uo)
+          shape_element_uo->executeJacobianWrapper(jvar, dof_indices);
+      }
+    }
+  }
+
   _fe_problem.swapBackMaterials(_tid);
 }
 

framework/src/base/FEProblem.C

@@ -2278,6 +2278,9 @@ FEProblem::execute(const ExecFlagType & exec_type)
 {
   // Set the current flag
   _current_execute_on_flag = exec_type;
+  if (exec_type == EXEC_NONLINEAR)
+    _currently_computing_jacobian = true;
+
 
   // Pre-aux UserObjects
   Moose::perf_log.push("computeUserObjects()", "Execution");
@@ -2301,6 +2304,7 @@ FEProblem::execute(const ExecFlagType & exec_type)
 
   // Return the current flag to None
   _current_execute_on_flag = EXEC_NONE;
+  _currently_computing_jacobian = false;
 }
 
 void

framework/src/userobject/ShapeElementUserObject.C

@@ -0,0 +1,58 @@
+/****************************************************************/
+/*               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 "ShapeElementUserObject.h"
+#include "Assembly.h"
+
+template<>
+InputParameters validParams<ShapeElementUserObject>()
+{
+  InputParameters params = validParams<ElementUserObject>();
+  params.addParam<bool>("compute_jacobians", true, "Compute Jacobians for coupled variables");
+  params.addParamNamesToGroup("compute_jacobians", "Advanced");
+  return params;
+}
+
+ShapeElementUserObject::ShapeElementUserObject(const InputParameters & parameters) :
+    ElementUserObject(parameters),
+    _phi(_assembly.phi()),
+    _grad_phi(_assembly.gradPhi()),
+    _currently_computing_jacobian(_fe_problem.currentlyComputingJacobian()),
+    _compute_jacobians(getParam<bool>("compute_jacobians"))
+{
+  mooseWarning("Jacobian calculation in UserObjects is an experimental capability with a potentially unstable interface.");
+}
+
+unsigned int
+ShapeElementUserObject::coupled(const std::string & var_name, unsigned int comp)
+{
+  MooseVariable * var = getVar(var_name, comp);
+
+  // add to the set of variables for which executeJacobian will be called
+  if (_compute_jacobians && var->kind() == Moose::VAR_NONLINEAR)
+    _jacobian_moose_variables.insert(var);
+
+  // return the variable number
+  return ElementUserObject::coupled(var_name, comp);
+}
+
+void
+ShapeElementUserObject::executeJacobianWrapper(unsigned int jvar, const std::vector<dof_id_type> & dof_indices)
+{
+  for (_j = 0; _j < _phi.size(); ++_j)
+  {
+    _j_global = dof_indices[_j];
+    executeJacobian(jvar);
+  }
+}

test/include/kernels/ExampleShapeElementKernel.h

@@ -0,0 +1,42 @@
+/****************************************************************/
+/*               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 EXAMPLESHAPEELEMENTKERNEL_H
+#define EXAMPLESHAPEELEMENTKERNEL_H
+
+#include "Kernel.h"
+#include "ExampleShapeElementUserObject.h"
+
+class ExampleShapeElementKernel : public Kernel
+{
+public:
+  ExampleShapeElementKernel(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpResidual();
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar);
+
+  const ExampleShapeElementUserObject & _shp;
+  const Real & _shp_integral;
+  const std::vector<Real> & _shp_jacobian;
+
+  unsigned int _u_var;
+  const std::vector<dof_id_type> & _u_dofs;
+  unsigned int _v_var;
+  const std::vector<dof_id_type> & _v_dofs;
+};
+
+template<>
+InputParameters validParams<ExampleShapeElementKernel>();
+
+#endif //EXAMPLESHAPEELEMENTKERNEL_H

test/include/userobjects/ExampleShapeElementUserObject.h

@@ -0,0 +1,64 @@
+/****************************************************************/
+/*               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 EXAMPLESHAPEELEMENTUSEROBJECT_H
+#define EXAMPLESHAPEELEMENTUSEROBJECT_H
+
+#include "ShapeElementUserObject.h"
+
+//Forward Declarations
+class ExampleShapeElementUserObject;
+
+template<>
+InputParameters validParams<ExampleShapeElementUserObject>();
+
+/**
+ * Test and proof of concept class for computing UserObject Jacobians using the
+ * ShapeElementUserObject base class. This object computes the integral
+ * \f$ \int_\Omega u^2v dr \f$
+ * and builds a vector of all derivatives of the integral w.r.t. the DOFs of u and v.
+ * These Jacobian terms can be utilized by a Kernel that uses the integral in the
+ * calculation of its residual.
+ */
+class ExampleShapeElementUserObject :
+  public ShapeElementUserObject
+{
+public:
+  ExampleShapeElementUserObject(const InputParameters & parameters);
+
+  virtual ~ExampleShapeElementUserObject() {}
+
+  virtual void initialize();
+  virtual void execute();
+  virtual void executeJacobian(unsigned int jvar);
+  virtual void finalize();
+  virtual void threadJoin(const UserObject & y);
+
+  ///@{ custom UserObject interface functions
+  const Real & getIntegral() const { return _integral; }
+  const std::vector<Real> & getJacobian() const { return _jacobian_storage; }
+  ///@}
+
+protected:
+  Real _integral;
+  std::vector<Real> _jacobian_storage;
+
+  const VariableValue & _u_value;
+  unsigned int _u_var;
+  const VariableValue & _v_value;
+  unsigned int _v_var;
+};
+
+#endif