Merge pull request #65 from aurora-multiphysics/alexanderianblair/h1-…

…integrators

Extend range of integrators and preconditioner types for use with H1 conforming FEs

doc/content/source/bcs/MFEMScalarBoundaryIntegratedBC.md

@@ -0,0 +1,21 @@
+# MFEMScalarBoundaryIntegratedBC
+
+## Summary
+
+!syntax description /BCs/MFEMScalarBoundaryIntegratedBC
+
+## Overview
+
+Adds the boundary integrator for integrating the linear form
+
+!equation
+(f, v)_{\partial\Omega} \,\,\, \forall v \in V
+
+where the test variable $v \in H^1$ and $f$ is a scalar coefficient. Often used for representing
+Neumann-type boundary conditions.
+
+!syntax parameters /BCs/MFEMScalarBoundaryIntegratedBC
+
+!syntax inputs /BCs/MFEMScalarBoundaryIntegratedBC
+
+!syntax children /BCs/MFEMScalarBoundaryIntegratedBC

doc/content/source/kernels/MFEMVectorDomainLFKernel.md

@@ -0,0 +1,28 @@
+# MFEMVectorDomainLFKernel
+
+!syntax description /Kernels/MFEMVectorDomainLFKernel
+
+## Overview
+
+Adds the domain integrator for integrating the linear form
+
+!equation
+(\vec f, \vec v)_\Omega \,\,\, \forall \vec v \in V
+
+where $\vec v \in \vec H^1$ is the test variable and $\vec f$ is a
+vector forcing coefficient.
+
+This term arises from the weak form of the forcing term
+
+!equation
+\vec f
+
+## Example Input File Syntax
+
+!listing kernels/gravity.i
+
+!syntax parameters /Kernels/MFEMVectorDomainLFKernel
+
+!syntax inputs /Kernels/MFEMVectorDomainLFKernel
+
+!syntax children /Kernels/MFEMVectorDomainLFKernel

doc/content/source/kernels/MFEMVectorFEDomainLFKernel.md

@@ -9,7 +9,7 @@ Adds the domain integrator for integrating the linear form
 !equation
 (\vec f, \vec v)_\Omega \,\,\, \forall \vec v \in V
 
-where $\vec v \in H(\mathrm{curl})$ or $H(\mathrm{div})$ is the trial variable and $\vec f$ is a
+where $\vec v \in H(\mathrm{curl})$ or $H(\mathrm{div})$ is the test variable and $\vec f$ is a
 vector forcing coefficient.
 
 This term arises from the weak form of the forcing term

include/bcs/MFEMScalarBoundaryIntegratedBC.h

@@ -0,0 +1,21 @@
+#pragma once
+#include "MFEMIntegratedBC.h"
+
+class MFEMScalarBoundaryIntegratedBC : public MFEMIntegratedBC
+{
+public:
+  static InputParameters validParams();
+
+  MFEMScalarBoundaryIntegratedBC(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 _coef_name;
+  mfem::Coefficient & _coef;
+};

include/kernels/MFEMVectorDomainLFKernel.h

@@ -0,0 +1,20 @@
+#pragma once
+#include "MFEMKernel.h"
+
+/*
+(\\vec f, \\vec u')
+*/
+class MFEMVectorDomainLFKernel : public MFEMKernel<mfem::LinearFormIntegrator>
+{
+public:
+  static InputParameters validParams();
+
+  MFEMVectorDomainLFKernel(const InputParameters & parameters);
+  ~MFEMVectorDomainLFKernel() override {}
+
+  virtual mfem::LinearFormIntegrator * createIntegrator() override;
+
+protected:
+  std::string _vec_coef_name;
+  mfem::VectorCoefficient & _vec_coef;
+};

include/solvers/MFEMHypreBoomerAMG.h

@@ -20,5 +20,7 @@ class MFEMHypreBoomerAMG : public MFEMSolverBase
   void constructSolver(const InputParameters & parameters) override;
 
 private:
+  std::shared_ptr<mfem::ParFiniteElementSpace> _mfem_fespace{nullptr};
+  mfem::real_t _strength_threshold;
   std::shared_ptr<mfem::HypreBoomerAMG> _solver{nullptr};
 };

src/bcs/MFEMScalarBoundaryIntegratedBC.C

@@ -0,0 +1,37 @@
+#include "MFEMScalarBoundaryIntegratedBC.h"
+
+registerMooseObject("PlatypusApp", MFEMScalarBoundaryIntegratedBC);
+
+InputParameters
+MFEMScalarBoundaryIntegratedBC::validParams()
+{
+  InputParameters params = MFEMIntegratedBC::validParams();
+  params.addClassDescription("Adds the domain integrator to an MFEM problem for the linear form "
+                             "$(f, v)_\\Omega$ "
+                             "arising from the weak form of the forcing term $f$.");
+  params.addRequiredParam<std::string>(
+      "coefficient", "The scalar MFEM coefficient which will be used in the integrated BC.");
+  return params;
+}
+
+MFEMScalarBoundaryIntegratedBC::MFEMScalarBoundaryIntegratedBC(const InputParameters & parameters)
+  : MFEMIntegratedBC(parameters),
+    _coef_name(getParam<std::string>("coefficient")),
+    _coef(getMFEMProblem().getProperties().getScalarProperty(_coef_name))
+{
+}
+
+// Create a new MFEM integrator to apply to the RHS of the weak form. Ownership managed by the
+// caller.
+mfem::LinearFormIntegrator *
+MFEMScalarBoundaryIntegratedBC::createLinearFormIntegrator()
+{
+  return new mfem::BoundaryLFIntegrator(_coef);
+}
+
+// Create a new MFEM integrator to apply to LHS of the weak form. Ownership managed by the caller.
+mfem::BilinearFormIntegrator *
+MFEMScalarBoundaryIntegratedBC::createBilinearFormIntegrator()
+{
+  return nullptr;
+}

src/kernels/MFEMVectorDomainLFKernel.C

@@ -0,0 +1,27 @@
+#include "MFEMVectorDomainLFKernel.h"
+
+registerMooseObject("PlatypusApp", MFEMVectorDomainLFKernel);
+
+InputParameters
+MFEMVectorDomainLFKernel::validParams()
+{
+  InputParameters params = MFEMKernel::validParams();
+  params.addClassDescription("Adds the domain integrator to an MFEM problem for the linear form "
+                             "$(\\vec f, \\vec v)_\\Omega$ "
+                             "arising from the weak form of the forcing term $\\vec f$.");
+  params.addParam<std::string>("vector_coefficient", "Name of body force density $\\vec f$.");
+  return params;
+}
+
+MFEMVectorDomainLFKernel::MFEMVectorDomainLFKernel(const InputParameters & parameters)
+  : MFEMKernel(parameters),
+    _vec_coef_name(getParam<std::string>("vector_coefficient")),
+    _vec_coef(getMFEMProblem().getProperties().getVectorProperty(_vec_coef_name))
+{
+}
+
+mfem::LinearFormIntegrator *
+MFEMVectorDomainLFKernel::createIntegrator()
+{
+  return new mfem::VectorDomainLFIntegrator(_vec_coef);
+}

src/solvers/MFEMHypreBoomerAMG.C

@@ -13,11 +13,18 @@ MFEMHypreBoomerAMG::validParams()
   params.addParam<double>("l_tol", 1e-5, "Set the relative tolerance.");
   params.addParam<int>("l_max_its", 10000, "Set the maximum number of iterations.");
   params.addParam<int>("print_level", 2, "Set the solver verbosity.");
+  params.addParam<UserObjectName>(
+      "fespace", "H1 FESpace to use in HypreBoomerAMG setup for elasticity problems.");
+  params.addParam<mfem::real_t>(
+      "strength_threshold", 0.25, "HypreBoomerAMG strong threshold. Defaults to 0.25.");
   return params;
 }
 
 MFEMHypreBoomerAMG::MFEMHypreBoomerAMG(const InputParameters & parameters)
-  : MFEMSolverBase(parameters)
+  : MFEMSolverBase(parameters),
+    _mfem_fespace(isParamSetByUser("fespace") ? getUserObject<MFEMFESpace>("fespace").getFESpace()
+                                              : nullptr),
+    _strength_threshold(getParam<mfem::real_t>("strength_threshold"))
 {
   constructSolver(parameters);
 }
@@ -30,4 +37,10 @@ MFEMHypreBoomerAMG::constructSolver(const InputParameters & parameters)
   _solver->SetTol(getParam<double>("l_tol"));
   _solver->SetMaxIter(getParam<int>("l_max_its"));
   _solver->SetPrintLevel(getParam<int>("print_level"));
+  _solver->SetStrengthThresh(_strength_threshold);
+
+  if (_mfem_fespace)
+  {
+    _solver->SetElasticityOptions(_mfem_fespace.get());
+  }
 }

test/tests/kernels/gold/OutputData/Gravity/Run0/Cycle000000/data.pvtu

@@ -0,0 +1,20 @@
+<?xml version="1.0"?>
+<VTKFile type="PUnstructuredGrid" version ="0.1" byte_order="LittleEndian">
+<PUnstructuredGrid GhostLevel="0">
+<PPoints>
+	<PDataArray type="Float64"  Name="Points" NumberOfComponents="3" format="ascii"/>
+</PPoints>
+<PCells>
+	<PDataArray type="Int32"  Name="connectivity" NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="Int32"  Name="offsets"      NumberOfComponents="1" format="ascii"/>
+	<PDataArray type="UInt8"  Name="types"        NumberOfComponents="1" format="ascii"/>
+</PCells>
+<PPointData>
+<PDataArray type="Float64" Name="displacement" NumberOfComponents="3" ComponentName0="0" ComponentName1="1" ComponentName2="2" format="ascii" />
+</PPointData>
+<PCellData>
+	<PDataArray type="Int32" Name="attribute" NumberOfComponents="1" format="ascii"/>
+</PCellData>
+<Piece Source="proc000000.vtu"/>
+</PUnstructuredGrid>
+</VTKFile>