Merge branch 'master' into integrate-nontensor-qf

fem/qspace.cpp

@@ -315,13 +315,36 @@ void FaceQuadratureSpace::Save(std::ostream &os) const
 
 const Vector &FaceQuadratureSpace::GetGeometricFactorWeights() const
 {
-   auto flags = FaceGeometricFactors::DETERMINANTS;
-   // TODO: assumes only one integration rule. This should be fixed once
-   // Mesh::GetFaceGeometricFactors acceps a QuadratureSpace instead of
-   // IntegrationRule.
-   const IntegrationRule &ir = GetIntRule(0);
-   auto *geom = mesh.GetFaceGeometricFactors(ir, flags, face_type);
-   return geom->detJ;
+   if (mesh.MeshGenerator() == (1 << 1)) // only tensor product elements
+   {
+      auto flags = FaceGeometricFactors::DETERMINANTS;
+      const IntegrationRule &ir = GetIntRule(0);
+      auto *geom = mesh.GetFaceGeometricFactors(ir, flags, face_type);
+      return geom->detJ;
+   }
+   else
+   {
+      if (non_tensor_weights.Size() == 0)
+      {
+         non_tensor_weights.SetSize(size);
+         int idx = 0;
+         for (int f = 0; f < face_indices.Size(); ++f)
+         {
+            const int f_idx = face_indices[f];
+            const Geometry::Type geom = mesh.GetFaceGeometry(f_idx);
+            ElementTransformation &T = *mesh.GetFaceTransformation(f_idx);
+            const IntegrationRule &ir = *int_rule[geom];
+            for (int q = 0; q < ir.Size(); ++q)
+            {
+               const IntegrationPoint &ip = ir.IntPoint(q);
+               T.SetIntPoint(&ip);
+               non_tensor_weights[idx] = T.Weight();
+               idx += 1;
+            }
+         }
+      }
+      return non_tensor_weights;
+   }
 }
 
 } // namespace mfem

fem/qspace.hpp

@@ -220,6 +220,9 @@ class FaceQuadratureSpace : public QuadratureSpaceBase
    /// Inverse of the map @a face_indices.
    const std::unordered_map<int,int> &face_indices_inv;
 
+   /// Integration weights (used only in the non-tensor case).
+   mutable Vector non_tensor_weights;
+
    const Vector &GetGeometricFactorWeights() const override;
    void ConstructOffsets();
    void Construct();

tests/unit/fem/test_quadf_coef.cpp

@@ -160,8 +160,10 @@ TEST_CASE("Quadrature Function Integration", "[QuadratureFunction][GPU]")
    auto fname = GENERATE(
                    "../../data/star.mesh",
                    "../../data/star-q3.mesh",
+                   "../../data/beam-tri.mesh",
                    "../../data/fichera.mesh",
-                   "../../data/fichera-q3.mesh"
+                   "../../data/fichera-q3.mesh",
+                   "../../data/beam-tet.mesh"
                 );
    const int order = GENERATE(1, 2, 3);