triangle test

examples/CMakeLists.txt

@@ -46,6 +46,7 @@ list(APPEND ALL_EXE_SRCS
   ex38.cpp
   ex39.cpp
   ex40.cpp
+  hmm.cpp
   )
 
 if (MFEM_USE_MPI)

examples/hmm.cpp

@@ -0,0 +1,139 @@
+// Hybrid mixed method
+
+#include "mfem.hpp"
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../data/inline-tri.mesh";
+   int order = 1;
+   int ref = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.AddOption(&ref, "-r", "--refine", "Refinement levels");
+   args.ParseCheck();
+
+   Mesh mesh(mesh_file);
+   for (int i = 0; i < ref; ++i) { mesh.UniformRefinement(); }
+
+   BrokenHdivFECollection fec_broken_rt(order - 1, mesh.Dimension());
+   L2_FECollection fec_l2(order - 1, mesh.Dimension());
+   DG_Interface_FECollection fec_trace(order - 1, mesh.Dimension());
+
+   FiniteElementSpace S_h(&mesh, &fec_broken_rt);
+   FiniteElementSpace V_h(&mesh, &fec_l2);
+   FiniteElementSpace M_h(&mesh, &fec_trace);
+
+   cout << "S_h size: " << S_h.GetTrueVSize() << '\n';
+   cout << "V_h size: " << V_h.GetTrueVSize() << '\n';
+   cout << "M_h size: " << M_h.GetTrueVSize() << '\n';
+
+   BilinearForm mass(&S_h);
+   mass.AddDomainIntegrator(new VectorFEMassIntegrator);
+   mass.Assemble();
+   mass.Finalize();
+
+   MixedBilinearForm div(&S_h, &V_h);
+   div.AddDomainIntegrator(new VectorFEDivergenceIntegrator);
+   div.Assemble();
+   div.Finalize();
+
+   MixedBilinearForm trace(&M_h, &S_h);
+   trace.AddTraceFaceIntegrator(new NormalTraceJumpIntegrator);
+   trace.Assemble();
+   trace.Finalize();
+
+   Array<int> empty;
+   Array<int> ess_trace_dofs;
+   M_h.GetBoundaryTrueDofs(ess_trace_dofs);
+
+   std::cout << "Boundary DOFs: " << ess_trace_dofs.Size() << '\n';
+
+   SparseMatrix T;
+   trace.FormRectangularSystemMatrix(ess_trace_dofs, empty, T);
+   unique_ptr<SparseMatrix> Tt(Transpose(T));
+
+   SparseMatrix &D(div.SpMat());
+   unique_ptr<SparseMatrix> Dt(Transpose(D));
+   (*Dt) *= -1.0;
+
+   SparseMatrix BC(M_h.GetTrueVSize(), M_h.GetTrueVSize());
+   for (int i : ess_trace_dofs)
+   {
+      BC.Set(i, i, 1.0);
+   }
+   BC.Finalize();
+
+   Array<int> offsets(4);
+   offsets[0] = 0;
+   offsets[1] = S_h.GetTrueVSize();
+   offsets[2] = offsets[1] + V_h.GetTrueVSize();
+   offsets[3] = offsets[2] + M_h.GetTrueVSize();
+
+   BlockMatrix matrix(offsets);
+   // Row 0
+   matrix.SetBlock(0, 0, &mass.SpMat());
+   matrix.SetBlock(0, 1, Dt.get());
+   matrix.SetBlock(0, 2, &T);
+   // Row 1
+   matrix.SetBlock(1, 0, &D);
+   // Row 2
+   matrix.SetBlock(2, 0, Tt.get());
+   matrix.SetBlock(2, 2, &BC);
+
+   unique_ptr<SparseMatrix> monolithic(matrix.CreateMonolithic());
+
+   // Form right-hand side
+   auto f = [](const Vector &xvec)
+   {
+      return 2*M_PI*M_PI*(sin(M_PI*xvec[0]) * sin(M_PI*xvec[1]));
+   };
+   FunctionCoefficient f_coeff(f);
+   LinearForm F(&V_h);
+   F.AddDomainIntegrator(new DomainLFIntegrator(f_coeff));
+   F.Assemble();
+
+   BlockVector rhs(offsets);
+   rhs.GetBlock(0) = 0.0;
+   rhs.GetBlock(1) = F;
+   rhs.GetBlock(2) = 0.0;
+
+   BlockVector solution(offsets);
+
+   UMFPackSolver solver(*monolithic);
+   solver.Mult(rhs, solution);
+
+   GridFunction sigma(&S_h);
+   GridFunction u(&V_h);
+   GridFunction lambda(&M_h);
+
+   sigma = solution.GetBlock(0);
+   u = solution.GetBlock(1);
+   lambda = solution.GetBlock(2);
+
+   auto u_exact = [](const Vector &xvec)
+   {
+      return sin(M_PI*xvec[0]) * sin(M_PI*xvec[1]);
+   };
+   FunctionCoefficient u_coeff(u_exact);
+
+   std::cout << "Error: " << u.ComputeL2Error(u_coeff) << '\n';
+
+   ParaViewDataCollection pv("HMM", &mesh);
+   pv.SetPrefixPath("ParaView");
+   pv.SetHighOrderOutput(true);
+   pv.SetLevelsOfDetail(order);
+   pv.RegisterField("u", &u);
+   pv.RegisterField("sigma", &sigma);
+   pv.SetCycle(0);
+   pv.SetTime(0.0);
+   pv.Save();
+
+   return 0;
+}

fem/fe_coll.cpp

@@ -2682,6 +2682,126 @@ RT_FECollection::~RT_FECollection()
    }
 }
 
+BrokenHdivFECollection::BrokenHdivFECollection(
+   const int order, const int dim, const int cb_type, const int ob_type)
+   : FiniteElementCollection(order + 1)
+   , dim(dim)
+   , cb_type(cb_type)
+   , ob_type(ob_type)
+{
+   int p = order;
+   MFEM_VERIFY(p >= 0, "BrokenHdivFECollection requires order >= 0.");
+
+   int cp_type = BasisType::GetQuadrature1D(cb_type);
+   int op_type = BasisType::GetQuadrature1D(ob_type);
+
+   if (Quadrature1D::CheckClosed(cp_type) == Quadrature1D::Invalid)
+   {
+      const char *cb_name = BasisType::Name(cb_type); // this may abort
+      MFEM_ABORT("unknown closed BasisType: " << cb_name);
+   }
+   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid &&
+       ob_type != BasisType::IntegratedGLL)
+   {
+      const char *ob_name = BasisType::Name(ob_type); // this may abort
+      MFEM_ABORT("unknown open BasisType: " << ob_name);
+   }
+
+   if (cb_type == BasisType::GaussLobatto &&
+       ob_type == BasisType::GaussLegendre)
+   {
+      snprintf(fec_name, 32, "BRT_%dD_P%d", dim, p);
+   }
+   else
+   {
+      snprintf(fec_name, 32, "BRT@%c%c_%dD_P%d",
+               (int)BasisType::GetChar(cb_type),
+               (int)BasisType::GetChar(ob_type), dim, p);
+   }
+
+   for (int g = 0; g < Geometry::NumGeom; g++)
+   {
+      RT_Elements[g] = nullptr;
+      RT_dof[g] = 0;
+   }
+
+   const int pp1 = p + 1;
+   const int pp2 = p + 2;
+   if (dim == 2)
+   {
+      // TODO: cb_type, ob_type for triangles
+      RT_Elements[Geometry::TRIANGLE] = new RT_TriangleElement(p);
+      RT_dof[Geometry::TRIANGLE] = pp1*(p+3);
+
+      RT_Elements[Geometry::SQUARE] = new RT_QuadrilateralElement(p, cb_type,
+                                                                  ob_type);
+      // two vector components * n_unk_face *
+      RT_dof[Geometry::SQUARE] = 2*pp1*pp2;
+   }
+   else if (dim == 3)
+   {
+      // TODO: cb_type, ob_type for tets
+      RT_Elements[Geometry::TETRAHEDRON] = new RT_TetrahedronElement(p);
+      RT_dof[Geometry::TETRAHEDRON] = pp1*pp2*(pp1 + 2)/2;
+
+      RT_Elements[Geometry::CUBE] = new RT_HexahedronElement(p, cb_type, ob_type);
+      RT_dof[Geometry::CUBE] = 3*pp1*pp2*pp2;
+
+      RT_Elements[Geometry::PRISM] = new RT_WedgeElement(p);
+      RT_dof[Geometry::PRISM] = pp1*pp2*(3*pp1 + 4)/2;
+
+      RT_Elements[Geometry::PYRAMID] = new RT0PyrFiniteElement(false);
+      RT_dof[Geometry::PYRAMID] = 0;
+   }
+   else
+   {
+      MFEM_ABORT("invalid dim = " << dim);
+   }
+}
+
+const FiniteElement *
+BrokenHdivFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
+   {
+      return RT_Elements[GeomType];
+   }
+   else
+   {
+      if (error_mode == RETURN_NULL) { return nullptr; }
+      MFEM_ABORT("RT Pyramid basis functions are not yet supported "
+                 "for order > 0.");
+      return NULL;
+   }
+}
+
+const int *BrokenHdivFECollection::DofOrderForOrientation(
+   Geometry::Type GeomType,
+   int Or) const
+{
+   if (GeomType == Geometry::SEGMENT)
+   {
+      return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
+   }
+   else if (GeomType == Geometry::TRIANGLE)
+   {
+      return TriDofOrd[Or%6];
+   }
+   else if (GeomType == Geometry::SQUARE)
+   {
+      return QuadDofOrd[Or%8];
+   }
+   return NULL;
+}
+
+BrokenHdivFECollection::~BrokenHdivFECollection()
+{
+   for (int g = 0; g < Geometry::NumGeom; g++)
+   {
+      delete RT_Elements[g];
+   }
+}
+
 
 RT_Trace_FECollection::RT_Trace_FECollection(const int p, const int dim,
                                              const int map_type,

fem/fe_coll.hpp

@@ -119,6 +119,8 @@ public:
    | ND@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(curl) | * | * / * | H_CURL | Nedelec vector elements |
    | ND_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0  | H_CURL | H^{1/2}-conforming trace elements for H(curl) defined on the interface between mesh elements (faces) |
    | ND_Trace@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0 | H_CURL | H^{1/2}-conforming trace elements for H(curl) defined on the interface between mesh elements (faces) |
+   | BRT_[DIM]_[ORDER] | L2 | * | 1 / 0 | H_DIV | Broken Raviart-Thomas vector elements |
+   | BRT@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | L2 | * | * / * | H_DIV | Broken Raviart-Thomas vector elements |
    | RT_[DIM]_[ORDER] | H(div) | * | 1 / 0 | H_DIV | Raviart-Thomas vector elements |
    | RT@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(div) | * | * / * | H_DIV | Raviart-Thomas vector elements |
    | RT_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0 | INTEGRAL | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
@@ -389,6 +391,45 @@ public:
    virtual ~L2_FECollection();
 };
 
+/// Broken H(div) collection (RT elements, no continuity)
+class BrokenHdivFECollection : public FiniteElementCollection
+{
+private:
+   int dim;
+   int cb_type; // closed BasisType
+   int ob_type; // open BasisType
+   char fec_name[32];
+   FiniteElement *RT_Elements[Geometry::NumGeom];
+   int RT_dof[Geometry::NumGeom];
+   int *SegDofOrd[2], *TriDofOrd[6], *QuadDofOrd[8];
+
+public:
+   BrokenHdivFECollection(const int p, const int dim,
+                          const int cb_type = BasisType::GaussLobatto,
+                          const int ob_type = BasisType::GaussLegendre);
+
+   const FiniteElement *
+   FiniteElementForGeometry(Geometry::Type GeomType) const override;
+
+   int DofForGeometry(Geometry::Type GeomType) const override
+   { return RT_dof[GeomType]; }
+
+   const int *DofOrderForOrientation(Geometry::Type GeomType,
+                                     int Or) const override;
+
+   const char *Name() const override { return fec_name; }
+
+   int GetContType() const override { return DISCONTINUOUS; }
+
+   int GetClosedBasisType() const { return cb_type; }
+   int GetOpenBasisType() const { return ob_type; }
+
+   FiniteElementCollection *Clone(int p) const override
+   { return new BrokenHdivFECollection(p, dim, cb_type, ob_type); }
+
+   virtual ~BrokenHdivFECollection();
+};
+
 /// Declare an alternative name for L2_FECollection = DG_FECollection
 typedef L2_FECollection DG_FECollection;