Some logic to be able to use BC info in GradFaceIntegrator (#198)

The idea here is to take advantage of all the data that has already
been set up in class BCIntegrator by making class GradFaceIntegrator
its friend.  Then we can use the internal std::unordered_map objects
in BCIntegrator, which map the boundary attribute to the appropriate
BoundaryCondition class, to get access to any required boundary
condition information.  In principle that is all that is necessary,
but the data we need it is not in a convenient form.  To address this,
the method computeBdrPrimitiveStateForGradient has been added to the
BoundaryCondition class.  In the base class, this method just sets the
boundary state to the internal state, which results in the same
(incorrect) behavior we have now.  In future commits, this method will
be implemented in the children of BoundaryCondition in order to
provide the right boundary state.

src/BCintegrator.hpp

@@ -47,6 +47,8 @@ using namespace mfem;
 
 // Boundary face term: <F.n(u),[w]>
 class BCintegrator : public NonlinearFormIntegrator {
+  friend class GradFaceIntegrator;
+
  protected:
   MPI_Groups *groupsMPI;
 

src/BoundaryCondition.cpp

@@ -52,6 +52,10 @@ BoundaryCondition::BoundaryCondition(RiemannSolver *_rsolver, GasMixture *_mixtu
 
 BoundaryCondition::~BoundaryCondition() {}
 
+void BoundaryCondition::computeBdrPrimitiveStateForGradient(const Vector &stateIn, Vector &stateBC) const {
+  stateBC = stateIn;
+}
+
 double BoundaryCondition::aggregateArea(int bndry_patchnum, MPI_Comm bc_comm) {
   double area = 0;
 

src/BoundaryCondition.hpp

@@ -76,6 +76,16 @@ class BoundaryCondition {
   virtual void computeBdrFlux(Vector &normal, Vector &stateIn, DenseMatrix &gradState, double radius,
                               Vector &bdrFlux) = 0;
 
+  /** \brief Set the boundary state used in the gradient evaluation
+   *
+   *  The jump in the state at the boundary appears in the
+   *  right-hand-side of the gradient solve.  This method sets this
+   *  boundary state to be equal to the interior state such that this
+   *  term is zero.  If that is not what you want, you must override
+   *  this method in a derived class.
+   */
+  virtual void computeBdrPrimitiveStateForGradient(const Vector &stateIn, Vector &stateBC) const;
+
   // holding function for any miscellaneous items needed to initialize BCs
   // prior to use (and require MPI)
   virtual void initBCs() = 0;

src/M2ulPhyS.cpp

@@ -541,17 +541,9 @@ void M2ulPhyS::initVariables() {
       cout << "Unknown ODE solver type: " << config.GetTimeIntegratorType() << '\n';
   }
 
-  //   gradUp_A = new ParNonlinearForm(gradUpfes);
-  //   gradUp_A->AddInteriorFaceIntegrator(
-  //       new GradFaceIntegrator(intRules, dim, num_equation) );
-  gradUp_A = new GradNonLinearForm(
-#ifdef _GPU_
-      vfes,
-#else
-      gradUpfes,
-#endif
-      intRules, dim, num_equation, gpuArrays, maxIntPoints, maxDofs);
+  gradUp_A = new GradNonLinearForm(gradUpfes, intRules, dim, num_equation, gpuArrays, maxIntPoints, maxDofs);
   gradUp_A->AddInteriorFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation));
+  gradUp_A->AddBdrFaceIntegrator(new GradFaceIntegrator(intRules, dim, num_equation, bcIntegrator));
 
   rhsOperator = new RHSoperator(iter, dim, num_equation, order, eqSystem, max_char_speed, intRules, intRuleType,
                                 fluxClass, mixture, d_mixture, chemistry_, transportPtr, radiation_, vfes, gpuArrays,

src/faceGradientIntegration.cpp

@@ -33,8 +33,9 @@
 #include "faceGradientIntegration.hpp"
 
 // Implementation of class FaceIntegrator
-GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation)
-    : dim(_dim), num_equation(_num_equation), intRules(_intRules) {}
+GradFaceIntegrator::GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation,
+                                       BCintegrator *bc)
+    : dim(_dim), num_equation(_num_equation), intRules(_intRules), bc_(bc) {}
 
 void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2,
                                             FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) {
@@ -76,18 +77,41 @@ void GradFaceIntegrator::AssembleFaceVector(const FiniteElement &el1, const Fini
   }
   const IntegrationRule *ir = &intRules->Get(Tr.GetGeometryType(), intorder);
 
+  // Quadrature point loop
   for (int i = 0; i < ir->GetNPoints(); i++) {
     const IntegrationPoint &ip = ir->IntPoint(i);
 
     Tr.SetAllIntPoints(&ip);  // set face and element int. points
 
     // Calculate basis functions on both elements at the face
     el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
-    el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+    if (Tr.Elem2No < 0) {
+      shape2 = shape1;
+    } else {
+      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+    }
 
     // Interpolate U values at the point
     elfun1_mat.MultTranspose(shape1, iUp1);
-    elfun2_mat.MultTranspose(shape2, iUp2);
+    if (Tr.Elem2No < 0) {
+      assert(bc_ != NULL);
+
+      const int attr = Tr.Attribute;
+      std::unordered_map<int, BoundaryCondition *>::const_iterator ibc = bc_->inletBCmap.find(attr);
+      std::unordered_map<int, BoundaryCondition *>::const_iterator obc = bc_->outletBCmap.find(attr);
+      std::unordered_map<int, BoundaryCondition *>::const_iterator wbc = bc_->wallBCmap.find(attr);
+      if (ibc != bc_->inletBCmap.end()) {
+        ibc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+      }
+      if (obc != bc_->outletBCmap.end()) {
+        obc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+      }
+      if (wbc != bc_->wallBCmap.end()) {
+        wbc->second->computeBdrPrimitiveStateForGradient(iUp1, iUp2);
+      }
+    } else {
+      elfun2_mat.MultTranspose(shape2, iUp2);
+    }
 
     // Compute average
     // NB: Code below is mathematically equivalent to mean = 0.5*(iUp1

src/faceGradientIntegration.hpp

@@ -34,6 +34,7 @@
 
 #include <tps_config.h>
 
+#include "BCintegrator.hpp"
 #include "tps_mfem_wrap.hpp"
 
 using namespace mfem;
@@ -44,9 +45,10 @@ class GradFaceIntegrator : public NonlinearFormIntegrator {
   const int dim;
   const int num_equation;
   IntegrationRules *intRules;
+  BCintegrator *bc_;  // NB: GradFaceIntegrator is a friend of BCintegrator
 
  public:
-  GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation);
+  GradFaceIntegrator(IntegrationRules *_intRules, const int _dim, const int _num_equation, BCintegrator *bc = NULL);
 
   virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr,
                                   const Vector &elfun, Vector &elvect);