Merge pull request #41 from aurora-multiphysics/k-collie/linear-elastic

Add linear elastic beam test and required functionality

.gitignore

@@ -223,6 +223,9 @@ config.status
 *MooseConfig.h.tmp
 conf_vars.mk
 
+# gdb history
+.gdb_history
+
 # GPerf performance
 *.prof
 

include/actions/SetMeshFESpaceAction.h

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "Action.h"
+#include "MFEMProblem.h"
+/**
+ * This class implements the action ensuring the mesh uses the same FE space as
+ * the displacement for mesh updates.
+ */
+class SetMeshFESpaceAction : public Action
+{
+public:
+  static InputParameters validParams();
+
+  SetMeshFESpaceAction(const InputParameters & parameters);
+
+  virtual void act() override;
+};

include/bcs/MFEMVectorBoundaryIntegratedBC.h

@@ -0,0 +1,22 @@
+#pragma once
+#include "MFEMIntegratedBC.h"
+#include "MFEMVectorFunctionCoefficient.h"
+
+class MFEMVectorBoundaryIntegratedBC : public MFEMIntegratedBC
+{
+public:
+  static InputParameters validParams();
+
+  MFEMVectorBoundaryIntegratedBC(const InputParameters & parameters);
+
+  // Create a new MFEM integrator to apply to the RHS of the weak form. Ownership managed by the
+  // caller.
+  virtual mfem::LinearFormIntegrator * createLinearFormIntegrator();
+
+  // Create a new MFEM integrator to apply to LHS of the weak form. Ownership managed by the caller.
+  virtual mfem::BilinearFormIntegrator * createBilinearFormIntegrator();
+
+protected:
+  std::string _vec_coef_name;
+  MFEMVectorCoefficient * _vec_coef;
+};

include/bcs/MFEMVectorDirichletBC.h

@@ -1,24 +1,11 @@
 #pragma once
-#include "MFEMEssentialBC.h"
-#include "MFEMVectorFunctionCoefficient.h"
 
-class MFEMVectorDirichletBC : public MFEMEssentialBC
-{
-protected:
-  MFEMVectorCoefficient * _vec_coef{nullptr};
-  enum APPLY_TYPE
-  {
-    STANDARD,
-    TANGENTIAL,
-    NORMAL
-  };
+#include "MFEMVectorDirichletBCBase.h"
 
+class MFEMVectorDirichletBC : public MFEMVectorDirichletBCBase
+{
 public:
-  static InputParameters validParams();
-
   MFEMVectorDirichletBC(const InputParameters & parameters);
-
+  ~MFEMVectorDirichletBC() override = default;
   void ApplyBC(mfem::GridFunction & gridfunc, mfem::Mesh * mesh_) override;
-
-  APPLY_TYPE _boundary_apply_type;
 };

include/bcs/MFEMVectorDirichletBCBase.h

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "MFEMEssentialBC.h"
+#include "MFEMVectorFunctionCoefficient.h"
+#include "boundary_conditions.h"
+
+class MFEMVectorDirichletBCBase : public MFEMEssentialBC
+{
+public:
+  static InputParameters validParams();
+
+  ~MFEMVectorDirichletBCBase() override = default;
+
+protected:
+  MFEMVectorDirichletBCBase(const InputParameters & parameters);
+  MFEMVectorCoefficient * _vec_coef{nullptr};
+};

include/bcs/MFEMVectorNormalIntegratedBC.h

@@ -18,5 +18,5 @@ class MFEMVectorNormalIntegratedBC : public MFEMIntegratedBC
 
 protected:
   std::string _vec_coef_name;
-  mfem::VectorCoefficient & _vec_coef;
+  MFEMVectorCoefficient * _vec_coef;
 };

include/bcs/MFEMVectorTangentialDirichletBC.h

@@ -0,0 +1,11 @@
+#pragma once
+
+#include "MFEMVectorDirichletBCBase.h"
+
+class MFEMVectorTangentialDirichletBC : public MFEMVectorDirichletBCBase
+{
+public:
+  MFEMVectorTangentialDirichletBC(const InputParameters & parameters);
+  ~MFEMVectorTangentialDirichletBC() override = default;
+  void ApplyBC(mfem::GridFunction & gridfunc, mfem::Mesh * mesh_) override;
+};

include/kernels/MFEMLinearElasticityKernel.h

@@ -0,0 +1,26 @@
+#pragma once
+#include "MFEMKernel.h"
+
+/*
+(cᵢₖⱼₗ∇uⱼ, ∇vᵢ),
+cᵢₖⱼₗ = λ δᵢₖδⱼₗ + μ (δᵢⱼδₖₗ + δᵢₗδⱼₖ),
+λ = (Eν)/((1-2ν)(1+ν)),
+μ = E/(2(1+ν)),
+E is Young's modulus,
+ν is Poisson's ratio
+*/
+class MFEMLinearElasticityKernel : public MFEMKernel<mfem::BilinearFormIntegrator>
+{
+public:
+  static InputParameters validParams();
+
+  MFEMLinearElasticityKernel(const InputParameters & parameters);
+
+  virtual mfem::BilinearFormIntegrator * createIntegrator() override;
+
+protected:
+  std::string _lambda_name;
+  std::string _mu_name;
+  mfem::Coefficient & _lambda;
+  mfem::Coefficient & _mu;
+};

include/mesh/MFEMMesh.h

@@ -49,6 +49,25 @@ class MFEMMesh : public FileMesh
    */
   std::unique_ptr<MooseMesh> safeClone() const override;
 
+  /**
+   * Returns true if mesh displacement is required.
+   */
+  bool shouldDisplace() const { return _mesh_displacement_variable.has_value(); }
+
+  /**
+   * Returns an optional reference to displacement variable name.
+   */
+  std::optional<std::reference_wrapper<std::string const>> getMeshDisplacementVariable() const
+  {
+    return _mesh_displacement_variable;
+  }
+
+  /**
+   * Displace the nodes of the mesh by the given displacement.
+   * Does not update FE spaces for variables.
+   */
+  void displace(mfem::GridFunction const & displacement);
+
 private:
   /**
    * Builds a placeholder mesh when no MOOSE mesh is required.
@@ -60,6 +79,11 @@ class MFEMMesh : public FileMesh
    */
   void uniformRefinement(mfem::Mesh & mesh, int nref);
 
+  /**
+   * Holds name of variable used for mesh displacement, if set.
+   */
+  std::optional<std::string> _mesh_displacement_variable;
+
   /**
    * Smart pointers to mfem::ParMesh object. Do not access directly.
    * Use the accessors instead.

include/problem/MFEMProblem.h

@@ -155,6 +155,17 @@ class MFEMProblem : public ExternalProblem
    */
   MFEMProblemData & getProblemData() { return _problem_data; }
 
+  /**
+   * Displace the mesh, if mesh displacement is enabled.
+   */
+  void displaceMesh();
+
+  /**
+   * Returns optional reference to the displacement GridFunction to apply to nodes.
+   */
+  std::optional<std::reference_wrapper<mfem::ParGridFunction const>>
+  getMeshDisplacementGridFunction();
+
 protected:
   /**
    * Template method for adding kernels. We can only add kernels using equation system problem

src/actions/SetMeshFESpaceAction.C

@@ -0,0 +1,36 @@
+#include "SetMeshFESpaceAction.h"
+
+registerMooseAction("PlatypusApp", SetMeshFESpaceAction, "set_mesh_fe_space");
+
+InputParameters
+SetMeshFESpaceAction::validParams()
+{
+  InputParameters params = Action::validParams();
+  params.addClassDescription("Set the mesh nodal finite element space to the same as the mesh "
+                             "displacement variable, if one is specified.");
+  return params;
+}
+
+SetMeshFESpaceAction::SetMeshFESpaceAction(const InputParameters & parameters) : Action(parameters)
+{
+}
+
+void
+SetMeshFESpaceAction::act()
+{
+  auto * mfem_problem = dynamic_cast<MFEMProblem *>(_problem.get());
+
+  if (!mfem_problem)
+  {
+    return;
+  }
+
+  mfem::ParMesh & mesh = mfem_problem->mesh().getMFEMParMesh();
+  auto const displacement = mfem_problem->getMeshDisplacementGridFunction();
+  if (!displacement)
+  {
+    return;
+  }
+
+  mesh.SetNodalFESpace(displacement.value().get().ParFESpace());
+}

src/base/PlatypusApp.C

@@ -1,3 +1,4 @@
+#include "ActionFactory.h"
 #include "PlatypusApp.h"
 #include "Moose.h"
 #include "AppFactory.h"
@@ -54,6 +55,11 @@ associateSyntaxInner(Syntax & syntax, ActionFactory & /*action_factory*/)
   addTaskDependency("add_elemental_field_variable", "add_mfem_fespaces");
   addTaskDependency("add_kernel", "add_mfem_fespaces");
 
+  // set mesh FE space
+  registerTask("set_mesh_fe_space", true);
+  addTaskDependency("set_mesh_fe_space", "add_variable");
+  addTaskDependency("set_mesh_fe_space", "init_mesh");
+
   // add preconditioning.
   registerMooseObjectTask("add_mfem_preconditioner", MFEMSolverBase, false);
   registerSyntaxTask("AddMFEMPreconditionerAction", "Preconditioner/*", "add_mfem_preconditioner");

src/bcs/MFEMVectorBoundaryIntegratedBC.C

@@ -0,0 +1,36 @@
+#include "MFEMVectorBoundaryIntegratedBC.h"
+
+registerMooseObject("PlatypusApp", MFEMVectorBoundaryIntegratedBC);
+
+InputParameters
+MFEMVectorBoundaryIntegratedBC::validParams()
+{
+  InputParameters params = MFEMIntegratedBC::validParams();
+  params.addRequiredParam<std::string>(
+      "vector_coefficient", "The vector MFEM coefficient which will be used in the integrated BC");
+  return params;
+}
+
+// TODO: Currently assumes the vector function coefficient is 3D
+MFEMVectorBoundaryIntegratedBC::MFEMVectorBoundaryIntegratedBC(const InputParameters & parameters)
+  : MFEMIntegratedBC(parameters),
+    _vec_coef_name(getParam<std::string>("vector_coefficient")),
+    _vec_coef(const_cast<MFEMVectorCoefficient *>(
+        &getUserObject<MFEMVectorCoefficient>("vector_coefficient")))
+{
+}
+
+// Create a new MFEM integrator to apply to the RHS of the weak form. Ownership managed by the
+// caller.
+mfem::LinearFormIntegrator *
+MFEMVectorBoundaryIntegratedBC::createLinearFormIntegrator()
+{
+  return new mfem::VectorBoundaryLFIntegrator(*_vec_coef->getVectorCoefficient());
+}
+
+// Create a new MFEM integrator to apply to LHS of the weak form. Ownership managed by the caller.
+mfem::BilinearFormIntegrator *
+MFEMVectorBoundaryIntegratedBC::createBilinearFormIntegrator()
+{
+  return nullptr;
+}

src/bcs/MFEMVectorDirichletBC.C

@@ -2,21 +2,8 @@
 
 registerMooseObject("PlatypusApp", MFEMVectorDirichletBC);
 
-InputParameters
-MFEMVectorDirichletBC::validParams()
-{
-  InputParameters params = MFEMEssentialBC::validParams();
-  params.addRequiredParam<UserObjectName>(
-      "vector_coefficient", "The vector MFEM coefficient to use in the Dirichlet condition");
-  return params;
-}
-
-// TODO: Currently assumes the vector function coefficient is 3D
 MFEMVectorDirichletBC::MFEMVectorDirichletBC(const InputParameters & parameters)
-  : MFEMEssentialBC(parameters),
-    _vec_coef(const_cast<MFEMVectorCoefficient *>(
-        &getUserObject<MFEMVectorCoefficient>("vector_coefficient"))),
-    _boundary_apply_type{APPLY_TYPE::TANGENTIAL}
+  : MFEMVectorDirichletBCBase(parameters)
 {
 }
 
@@ -25,21 +12,5 @@ MFEMVectorDirichletBC::ApplyBC(mfem::GridFunction & gridfunc, mfem::Mesh * mesh_
 {
   mfem::Array<int> ess_bdrs(mesh_->bdr_attributes.Max());
   ess_bdrs = GetMarkers(*mesh_);
-  if (_vec_coef == nullptr)
-  {
-    MFEM_ABORT("Boundary condition does not store valid coefficients to specify the "
-               "components of the vector at the Dirichlet boundary.");
-  }
-
-  switch (_boundary_apply_type)
-  {
-    case STANDARD:
-      gridfunc.ProjectBdrCoefficient(*_vec_coef->getVectorCoefficient(), ess_bdrs);
-      break;
-    case NORMAL:
-      gridfunc.ProjectBdrCoefficientNormal(*_vec_coef->getVectorCoefficient(), ess_bdrs);
-      break;
-    case TANGENTIAL:
-      gridfunc.ProjectBdrCoefficientTangent(*_vec_coef->getVectorCoefficient(), ess_bdrs);
-  }
-}
+  gridfunc.ProjectBdrCoefficient(*_vec_coef->getVectorCoefficient(), ess_bdrs);
+}

src/bcs/MFEMVectorDirichletBCBase.C

@@ -0,0 +1,18 @@
+#include "MFEMVectorDirichletBCBase.h"
+
+InputParameters
+MFEMVectorDirichletBCBase::validParams()
+{
+  InputParameters params = MFEMEssentialBC::validParams();
+  params.addRequiredParam<UserObjectName>(
+      "vector_coefficient", "The vector MFEM coefficient to use in the Dirichlet condition");
+  return params;
+}
+
+// TODO: Currently assumes the vector function coefficient is 3D
+MFEMVectorDirichletBCBase::MFEMVectorDirichletBCBase(const InputParameters & parameters)
+  : MFEMEssentialBC(parameters),
+    _vec_coef(const_cast<MFEMVectorCoefficient *>(
+        &getUserObject<MFEMVectorCoefficient>("vector_coefficient")))
+{
+}

src/bcs/MFEMVectorNormalIntegratedBC.C

@@ -16,7 +16,8 @@ MFEMVectorNormalIntegratedBC::validParams()
 MFEMVectorNormalIntegratedBC::MFEMVectorNormalIntegratedBC(const InputParameters & parameters)
   : MFEMIntegratedBC(parameters),
     _vec_coef_name(getParam<std::string>("vector_coefficient")),
-    _vec_coef(getMFEMProblem().getProperties().getVectorProperty(_vec_coef_name))
+    _vec_coef(const_cast<MFEMVectorCoefficient *>(
+        &getUserObject<MFEMVectorCoefficient>("vector_coefficient")))
 {
 }
 
@@ -25,12 +26,12 @@ MFEMVectorNormalIntegratedBC::MFEMVectorNormalIntegratedBC(const InputParameters
 mfem::LinearFormIntegrator *
 MFEMVectorNormalIntegratedBC::createLinearFormIntegrator()
 {
-  return new mfem::BoundaryNormalLFIntegrator(_vec_coef);
+  return new mfem::BoundaryNormalLFIntegrator(*_vec_coef->getVectorCoefficient());
 }
 
 // Create a new MFEM integrator to apply to LHS of the weak form. Ownership managed by the caller.
 mfem::BilinearFormIntegrator *
 MFEMVectorNormalIntegratedBC::createBilinearFormIntegrator()
 {
   return nullptr;
-}
+}