Merge branch 'master' into matrix-free-FCT

Conflicts:
	examples/ex9.cpp
	fem/fe.cpp

data/unit-square.mesh

@@ -0,0 +1,97 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+9
+1 3 0 1 5 4
+2 3 1 2 6 5
+3 3 2 3 7 6
+4 3 4 5 9 8
+5 3 5 6 10 9
+6 3 6 7 11 10
+7 3 8 9 13 12
+8 3 9 10 14 13
+9 3 10 11 15 14
+
+boundary
+12
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 7
+5 1 7 11
+6 1 11 15
+7 1 15 14
+8 1 14 13
+9 1 13 12
+10 1 12 8
+11 1 8 4
+12 1 4 0
+
+vertices
+16
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 2
+Ordering: 1
+
+0 0
+0.333333333 0
+0 0.333333333
+0.333333333 0.333333333
+
+0.333333333 0
+0.666666667 0
+0.333333333 0.333333333
+0.666666667 0.333333333
+
+0.666666667 0
+1 0
+0.666666667 0.333333333
+1 0.333333333
+
+0 0.333333333
+0.333333333 0.333333333
+0 0.666666667
+0.333333333 0.666666667
+
+0.333333333 0.333333333
+0.666666667 0.333333333
+0.333333333 0.666666667
+0.666666667 0.666666667
+
+0.666666667 0.333333333
+1 0.333333333
+0.666666667 0.666666667
+1 0.666666667
+
+0 0.666666667
+0.333333333 0.666666667
+0 1
+0.333333333 1
+
+0.333333333 0.666666667
+0.666666667 0.666666667
+0.333333333 1
+0.666666667 1
+
+0.666666667 0.666666667
+1 0.666666667
+0.666666667 1
+1 1

fem/bilininteg.cpp

@@ -900,6 +900,52 @@ void ConvectionIntegrator::AssembleElementMatrix(
 }
 
 
+void MixedConvectionIntegrator::AssembleElementMatrix2(
+   const FiniteElement &tr_el, const FiniteElement &te_el,
+   ElementTransformation &Trans, DenseMatrix &elmat)
+{
+   int tr_nd = tr_el.GetDof();
+   int te_nd = te_el.GetDof();
+   int dim = te_el.GetDim(); // Using test geometry.
+
+#ifdef MFEM_THREAD_SAFE
+   DenseMatrix dshape, adjJ, Q_ir;
+   Vector shape, vec2, BdFidxT;
+#endif
+   elmat.SetSize(te_nd, tr_nd);
+   dshape.SetSize(tr_nd,dim);
+   adjJ.SetSize(dim);
+   shape.SetSize(te_nd);
+   vec2.SetSize(dim);
+   BdFidxT.SetSize(tr_nd);
+
+   Vector vec1;
+
+   // Using midpoint rule and test geometry.
+   const IntegrationRule *ir = &IntRules.Get(te_el.GetGeomType(), 1);
+
+   Q.Eval(Q_ir, Trans, *ir);
+
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      tr_el.CalcDShape(ip, dshape);
+      te_el.CalcShape(ip, shape);
+
+      Trans.SetIntPoint(&ip);
+      CalcAdjugate(Trans.Jacobian(), adjJ);
+      Q_ir.GetColumnReference(i, vec1);
+      vec1 *= alpha * ip.weight;
+
+      adjJ.Mult(vec1, vec2);
+      dshape.Mult(vec2, BdFidxT);
+
+      AddMultVWt(shape, BdFidxT, elmat);
+   }
+}
+
+
 void GroupConvectionIntegrator::AssembleElementMatrix(
    const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
 {
@@ -3428,4 +3474,63 @@ VectorInnerProductInterpolator::AssembleElementMatrix2(
    ran_fe.Project(dom_shape_coeff, Trans, elmat_as_vec);
 }
 
+
+void PrecondConvectionIntegrator::AssembleElementMatrix(
+   const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
+{
+   int i, nd = el.GetDof(), dim = el.GetDim();
+
+#ifdef MFEM_THREAD_SAFE
+   DenseMatrix dshape, adjJ, Q_ir;
+   Vector shape, vec2, BdFidxT;
+#endif
+   elmat.SetSize(nd);
+   dshape.SetSize(nd,dim);
+   adjJ.SetSize(dim);
+   shape.SetSize(nd);
+   vec2.SetSize(dim);
+   BdFidxT.SetSize(nd);
+
+   double w;
+   Vector vec1;
+   DenseMatrix mass(nd,nd), conv(nd,nd), lumpedM(nd,nd), tmp(nd,nd);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = Trans.OrderGrad(&el) + Trans.Order() + el.GetOrder();
+      order = max(order, 2 * el.GetOrder() + Trans.OrderW());
+      ir = &IntRules.Get(el.GetGeomType(), order);
+   }
+
+   Q.Eval(Q_ir, Trans, *ir);
+
+   conv = mass = 0.0;
+   for (i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      el.CalcDShape(ip, dshape);
+      el.CalcShape(ip, shape);
+
+      Trans.SetIntPoint(&ip);
+      CalcAdjugate(Trans.Jacobian(), adjJ);
+      Q_ir.GetColumnReference(i, vec1);
+      vec1 *= alpha * ip.weight;
+
+      adjJ.Mult(vec1, vec2);
+      dshape.Mult(vec2, BdFidxT);
+
+      AddMultVWt(shape, BdFidxT, conv);
+
+      w = Trans.Weight() * ip.weight;
+      AddMult_a_VVt(w, shape, mass);
+   }
+   lumpedM = mass;
+   lumpedM.Lump();
+   mass.Invert();
+
+   MultABt(mass, lumpedM, tmp);
+   MultAtB(tmp, conv, elmat); // using symmetry of mass matrix
+}
+
 }

fem/bilininteg.hpp

@@ -1731,6 +1731,26 @@ public:
                                       DenseMatrix &);
 };
 
+/// alpha (q . grad u, v)
+class MixedConvectionIntegrator : public BilinearFormIntegrator
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   DenseMatrix dshape, adjJ, Q_ir;
+   Vector shape, vec2, BdFidxT;
+#endif
+   VectorCoefficient &Q;
+   double alpha;
+
+public:
+   MixedConvectionIntegrator(VectorCoefficient &q, double a = 1.0)
+      : Q(q) { alpha = a; }
+   virtual void AssembleElementMatrix2(const FiniteElement &,
+                                       const FiniteElement &,
+                                       ElementTransformation &,
+                                       DenseMatrix &);
+};
+
 /// alpha (q . grad u, v) using the "group" FE discretization
 class GroupConvectionIntegrator : public BilinearFormIntegrator
 {
@@ -2494,6 +2514,28 @@ protected:
    VectorCoefficient &VQ;
 };
 
+
+/** Class for local assembly of M_L M_C^-1 K, where M_L and M_C are
+    the lumped and consistent mass matrices and K is the convection
+    matrix. The spaces are assumed to be L2 conforming. */
+class PrecondConvectionIntegrator: public BilinearFormIntegrator
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   DenseMatrix dshape, adjJ, Q_ir;
+   Vector shape, vec2, BdFidxT;
+#endif
+   VectorCoefficient &Q;
+   double alpha;
+
+public:
+   PrecondConvectionIntegrator(VectorCoefficient &q, double a = 1.0)
+      : Q(q) { alpha = a; }
+   virtual void AssembleElementMatrix(const FiniteElement &,
+                                      ElementTransformation &,
+                                      DenseMatrix &);
+};
+
 }
 
 #endif

fem/fe.cpp

@@ -203,6 +203,12 @@ void FiniteElement::CalcPhysDShape(ElementTransformation &Trans,
    Mult(vshape, Trans.InverseJacobian(), dshape);
 }
 
+void FiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   FiniteElements!");
+}
+
 
 void ScalarFiniteElement::NodalLocalInterpolation (
    ElementTransformation &Trans, DenseMatrix &I,
@@ -278,6 +284,12 @@ void ScalarFiniteElement::ScalarLocalInterpolation(
    }
 }
 
+void ScalarFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   ScalarFiniteElements!");
+}
+
 
 void NodalFiniteElement::ProjectCurl_2D(
    const FiniteElement &fe, ElementTransformation &Trans,
@@ -505,6 +517,12 @@ void NodalFiniteElement::ProjectDiv(
    }
 }
 
+void NodalFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   NodalFiniteElements!");
+}
+
 
 void PositiveFiniteElement::Project(
    Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
@@ -561,6 +579,11 @@ void PositiveFiniteElement::Project(
    }
 }
 
+void PositiveFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   PositiveFiniteElements!");
+}
 
 void VectorFiniteElement::CalcShape (
    const IntegrationPoint &ip, Vector &shape ) const
@@ -1002,6 +1025,12 @@ void VectorFiniteElement::LocalInterpolation_ND(
    }
 }
 
+void VectorFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   VectorFiniteElements!");
+}
+
 void VectorFiniteElement::LocalRestriction_RT(
    const double *nk, const Array<int> &d2n, ElementTransformation &Trans,
    DenseMatrix &R) const
@@ -1386,6 +1415,14 @@ Quad2DFiniteElement::Quad2DFiniteElement()
    Nodes.IntPoint(5).y = 0.5;
 }
 
+void QuadPos1DFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   dofs.SetSize(1,2);
+   dofs(0,0) = 0;
+   dofs(0,1) = 1;
+}
+
+
 void Quad2DFiniteElement::CalcShape(const IntegrationPoint &ip,
                                     Vector &shape) const
 {
@@ -1819,6 +1856,15 @@ void BiQuadPos2DFiniteElement::Project (
    }
 }
 
+void BiQuadPos2DFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   dofs.SetSize(3,4);
+   dofs(0,0) = 0;  dofs(1,0) = 4;  dofs(2,0) = 1;
+   dofs(0,1) = 1;  dofs(1,1) = 5;  dofs(2,1) = 2;
+   dofs(0,2) = 2;  dofs(1,2) = 6;  dofs(2,2) = 3;
+   dofs(0,3) = 3;  dofs(1,3) = 7;  dofs(2,3) = 0;
+}
+
 
 GaussBiQuad2DFiniteElement::GaussBiQuad2DFiniteElement()
    : NodalFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
@@ -7071,6 +7117,12 @@ PositiveTensorFiniteElement::PositiveTensorFiniteElement(
                            dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
      TensorBasisElement(dims, p, BasisType::Positive, dmtype) { }
 
+void PositiveTensorFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   PositiveTensorFiniteElements!");
+}
+
 
 H1_SegmentElement::H1_SegmentElement(const int p, const int btype)
    : NodalTensorFiniteElement(1, p, VerifyClosed(btype), H1_DOF_MAP)
@@ -7489,6 +7541,13 @@ void H1Pos_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
    dofs[vertex] = 1.0;
 }
 
+void H1Pos_SegmentElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   dofs.SetSize(1,2);
+   dofs(0,0) = 0;
+   dofs(0,1) = 1;
+}
+
 
 H1Pos_QuadrilateralElement::H1Pos_QuadrilateralElement(const int p)
    : PositiveTensorFiniteElement(2, p, H1_DOF_MAP)
@@ -7557,6 +7616,19 @@ void H1Pos_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
    dofs[vertex] = 1.0;
 }
 
+void H1Pos_QuadrilateralElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int p = Order;
+   dofs.SetSize(p+1,4);
+   for (int i = 0; i <= p; i++)
+   {
+      dofs(i,0) = i;
+      dofs(i,1) = i*(p+1) + p;
+      dofs(i,2) = (p+1)*(p+1) - 1 - i;
+      dofs(i,3) = (p-i)*(p+1);
+   }
+}
+
 
 H1Pos_HexahedronElement::H1Pos_HexahedronElement(const int p)
    : PositiveTensorFiniteElement(3, p, H1_DOF_MAP)
@@ -7631,6 +7703,57 @@ void H1Pos_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
    dofs[vertex] = 1.0;
 }
 
+void H1Pos_HexahedronElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int p = Order;
+   dofs.SetSize((p+1)*(p+1), 6);
+   for (int bdrID = 0; bdrID < 6; bdrID++)
+   {
+      int o(0);
+      switch (bdrID)
+      {
+         case 0:
+            for (int i = 0; i < (p+1)*(p+1); i++)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+         case 1:
+            for (int i = 0; i <= p*(p+1)*(p+1); i+=(p+1)*(p+1))
+               for (int j = 0; j < p+1; j++)
+               {
+                  dofs(o++,bdrID) = i+j;
+               }
+            break;
+         case 2:
+            for (int i = p; i < (p+1)*(p+1)*(p+1); i+=p+1)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+         case 3:
+            for (int i = 0; i <= p*(p+1)*(p+1); i+=(p+1)*(p+1))
+               for (int j = p*(p+1); j < (p+1)*(p+1); j++)
+               {
+                  dofs(o++,bdrID) = i+j;
+               }
+            break;
+         case 4:
+            for (int i = 0; i <= (p+1)*((p+1)*(p+1)-1); i+=p+1)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+         case 5:
+            for (int i = p*(p+1)*(p+1); i < (p+1)*(p+1)*(p+1); i++)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+      }
+   }
+}
+
 
 H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
    : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
@@ -8151,6 +8274,19 @@ void H1Pos_TriangleElement::CalcDShape(const IntegrationPoint &ip,
    }
 }
 
+void H1Pos_TriangleElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int ctr = 0, p = Order;
+   dofs.SetSize(p+1, 3);
+   for (int i = 0; i <= p; i++)
+   {
+      dofs(i,0) = i;
+      dofs(i,1) = ctr + p;
+      dofs(i,2) = ctr + i;
+      ctr += p - i;
+   }
+}
+
 
 H1Pos_TetrahedronElement::H1Pos_TetrahedronElement(const int p)
    : PositiveFiniteElement(3, Geometry::TETRAHEDRON,
@@ -8403,6 +8539,59 @@ void H1Pos_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
    }
 }
 
+void H1Pos_TetrahedronElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int ctr, p = Order;
+   dofs.SetSize((p+1)*(p+2)/2, 4);
+   for (int bdrID = 0; bdrID < 4; bdrID++)
+   {
+      int o = 0;
+      switch (bdrID)
+      {
+         case 0:
+            ctr = p;
+            for (int i = 0; i <= p; i++)
+            {
+               for (int j = 0; j <= p - i; j++)
+               {
+                  dofs(o++,bdrID) = ctr;
+                  ctr += p - i - j;
+               }
+               ctr += p - i;
+            }
+            break;
+         case 1:
+            ctr = 0;
+            for (int i = 0; i <= p; i++)
+            {
+               for (int j = 0; j <= p - i; j++)
+               {
+                  dofs(o++,bdrID) = ctr;
+                  ctr += p + 1 - i - j;
+               }
+            }
+            break;
+         case 2:
+            ctr = 0;
+            for (int i = 0; i <= p; i++)
+            {
+               for (int j = 0; j <= p - i; j++)
+               {
+                  dofs(o++,bdrID) = ctr++;
+               }
+               ctr += - p + i - 1 + (p-i+1)*(p-i+2)/2;
+            }
+            break;
+         case 3:
+            for (int i = 0; i < (p+1)*(p+2)/2; i++)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+      }
+   }
+}
+
 
 H1_WedgeElement::H1_WedgeElement(const int p,
                                  const int btype)
@@ -8793,6 +8982,13 @@ void L2Pos_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
    dofs[vertex*Order] = 1.0;
 }
 
+void L2Pos_SegmentElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   dofs.SetSize(1,2);
+   dofs(0,0) = 0;
+   dofs(0,1) = Order;
+}
+
 
 L2_QuadrilateralElement::L2_QuadrilateralElement(const int p, const int btype)
    : NodalTensorFiniteElement(2, p, VerifyOpen(btype), L2_DOF_MAP)
@@ -8978,6 +9174,19 @@ void L2Pos_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
    }
 }
 
+void L2Pos_QuadrilateralElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int p = Order;
+   dofs.SetSize(p+1,4);
+   for (int i = 0; i <= p; i++)
+   {
+      dofs(i,0) = i;
+      dofs(i,1) = i*(p+1) + p;
+      dofs(i,2) = (p+1)*(p+1) - 1 - i;
+      dofs(i,3) = (p-i)*(p+1);
+   }
+}
+
 
 L2_HexahedronElement::L2_HexahedronElement(const int p, const int btype)
    : NodalTensorFiniteElement(3, p, VerifyOpen(btype), L2_DOF_MAP)
@@ -9221,6 +9430,57 @@ void L2Pos_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
    }
 }
 
+void L2Pos_HexahedronElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int p = Order;
+   dofs.SetSize((p+1)*(p+1), 6);
+   for (int bdrID = 0; bdrID < 6; bdrID++)
+   {
+      int o(0);
+      switch (bdrID)
+      {
+         case 0:
+            for (int i = 0; i < (p+1)*(p+1); i++)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+         case 1:
+            for (int i = 0; i <= p*(p+1)*(p+1); i+=(p+1)*(p+1))
+               for (int j = 0; j < p+1; j++)
+               {
+                  dofs(o++,bdrID) = i+j;
+               }
+            break;
+         case 2:
+            for (int i = p; i < (p+1)*(p+1)*(p+1); i+=p+1)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+         case 3:
+            for (int i = 0; i <= p*(p+1)*(p+1); i+=(p+1)*(p+1))
+               for (int j = p*(p+1); j < (p+1)*(p+1); j++)
+               {
+                  dofs(o++,bdrID) = i+j;
+               }
+            break;
+         case 4:
+            for (int i = 0; i <= (p+1)*((p+1)*(p+1)-1); i+=p+1)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+         case 5:
+            for (int i = p*(p+1)*(p+1); i < (p+1)*(p+1)*(p+1); i++)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+      }
+   }
+}
+
 
 L2_TriangleElement::L2_TriangleElement(const int p, const int btype)
    : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
@@ -9397,6 +9657,19 @@ void L2Pos_TriangleElement::ProjectDelta(int vertex, Vector &dofs) const
    }
 }
 
+void L2Pos_TriangleElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int ctr = 0, p = Order;
+   dofs.SetSize(p+1, 3);
+   for (int i = 0; i <= p; i++)
+   {
+      dofs(i,0) = i;
+      dofs(i,1) = ctr + p;
+      dofs(i,2) = ctr + i;
+      ctr += p - i;
+   }
+}
+
 
 L2_TetrahedronElement::L2_TetrahedronElement(const int p, const int btype)
    : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
@@ -9594,6 +9867,59 @@ void L2Pos_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
    }
 }
 
+void L2Pos_TetrahedronElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   int ctr, p = Order;
+   dofs.SetSize((p+1)*(p+2)/2, 4);
+   for (int bdrID = 0; bdrID < 4; bdrID++)
+   {
+      int o = 0;
+      switch (bdrID)
+      {
+         case 0:
+            ctr = p;
+            for (int i = 0; i <= p; i++)
+            {
+               for (int j = 0; j <= p - i; j++)
+               {
+                  dofs(o++,bdrID) = ctr;
+                  ctr += p - i - j;
+               }
+               ctr += p - i;
+            }
+            break;
+         case 1:
+            ctr = 0;
+            for (int i = 0; i <= p; i++)
+            {
+               for (int j = 0; j <= p - i; j++)
+               {
+                  dofs(o++,bdrID) = ctr;
+                  ctr += p + 1 - i - j;
+               }
+            }
+            break;
+         case 2:
+            ctr = 0;
+            for (int i = 0; i <= p; i++)
+            {
+               for (int j = 0; j <= p - i; j++)
+               {
+                  dofs(o++,bdrID) = ctr++;
+               }
+               ctr += - p + i - 1 + (p-i+1)*(p-i+2)/2;
+            }
+            break;
+         case 3:
+            for (int i = 0; i < (p+1)*(p+2)/2; i++)
+            {
+               dofs(o++,bdrID) = i;
+            }
+            break;
+      }
+   }
+}
+
 
 L2_WedgeElement::L2_WedgeElement(const int p, const int btype)
    : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
@@ -11644,6 +11970,14 @@ void ND_SegmentElement::CalcVShape(const IntegrationPoint &ip,
    obasis1d.Eval(ip.x, vshape);
 }
 
+
+void NURBSFiniteElement::ExtractBdrDofs(DenseMatrix &dofs) const
+{
+   mfem_error ("Error: Cannot use ExtractBdrDofs(...) function with\n"
+               "   NURBSFiniteElements!");
+}
+
+
 void NURBS1DFiniteElement::SetOrder() const
 {
    Order = kv[0]->GetOrder();

fem/fe.hpp

@@ -448,6 +448,12 @@ public:
    {
       return BasisType::CheckNodal(b_type);
    }
+
+   /** Routine that extracts the indices of all p-th order Bernstein basis functions
+       that are non-zero on each of the (dim-1)-dimensional boundaries of a finite
+       element. The columns of dofs hold the indices of basis functions corresponding
+       to one respective boundary defined according to the class Geometry. */
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 class ScalarFiniteElement : public FiniteElement
@@ -494,6 +500,9 @@ public:
    void ScalarLocalInterpolation(ElementTransformation &Trans,
                                  DenseMatrix &I,
                                  const ScalarFiniteElement &fine_fe) const;
+
+   /// Overrides the ExtractBdrDofs function to print an error.
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 class NodalFiniteElement : public ScalarFiniteElement
@@ -540,6 +549,9 @@ public:
    virtual void ProjectDiv(const FiniteElement &fe,
                            ElementTransformation &Trans,
                            DenseMatrix &div) const;
+
+   /// Overrides the ExtractBdrDofs function to print an error.
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -572,6 +584,9 @@ public:
 
    virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
                         DenseMatrix &I) const;
+
+   /// Overrides the ExtractBdrDofs function to print an error.
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 class VectorFiniteElement : public FiniteElement
@@ -679,6 +694,9 @@ public:
       FiniteElement(D, G, Do, O, F), Jinv(D)
    { RangeType = VECTOR; MapType = M; SetDerivMembers(); }
 #endif
+
+   /// Overrides the ExtractBdrDofs function to print an error.
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 class PointFiniteElement : public NodalFiniteElement
@@ -821,6 +839,7 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 /// Class for quadratic FE on triangle
@@ -900,6 +919,7 @@ public:
                         Vector &dofs) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const
    { dofs = 0.; dofs(vertex) = 1.; }
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 /// Bi-quadratic element on quad with nodes at the 9 Gaussian points
@@ -1758,6 +1778,9 @@ class PositiveTensorFiniteElement : public PositiveFiniteElement,
 public:
    PositiveTensorFiniteElement(const int dims, const int p,
                                const DofMapType dmtype);
+
+   /// Overrides the ExtractBdrDofs function to print an error.
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 class H1_SegmentElement : public NodalTensorFiniteElement
@@ -1826,6 +1849,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -1843,6 +1867,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -1860,6 +1885,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -1928,6 +1954,7 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -1954,6 +1981,7 @@ public:
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -2051,6 +2079,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -2088,6 +2117,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -2121,6 +2151,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -2160,6 +2191,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -2196,6 +2228,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 
@@ -2750,6 +2783,9 @@ public:
    Vector              &Weights    ()         const { return weights; }
    /// Update the NURBSFiniteElement according to the currently set knot vectors
    virtual void         SetOrder   ()         const { }
+
+   /// Overrides the ExtractBdrDofs function to print an error.
+   virtual void ExtractBdrDofs(DenseMatrix &dofs) const;
 };
 
 class NURBS1DFiniteElement : public NURBSFiniteElement

linalg/densemat.cpp

@@ -3917,7 +3917,6 @@ void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
    }
 }
 
-
 void LUFactors::Factor(int m)
 {
 #ifdef MFEM_USE_LAPACK