Implementation of MINI elements

[lmolin3]

Hi all,
I was trying to implement MINI elements, to solve a Navier-Stokes problem with P1bubble-P1 pair instead of Taylor-Hood P2-P1 elements.
I created a new FiniteElement adding the bubble function in the center for both Triangle and Tetrahedron geometries, and created the associated FiniteElementCollection class.
When running a simple test for Poiseulle flow on a rectangular domain I still get spurious pressure mode similar to the P1-P1 case, even when using P1bubble-P1, which are supposed to be stable.
Any idea what might be missing?

Thanks for any help!

Leo

I'm attaching the code for both the FiniteElement and FiniteElementCollection.

The FINITE ELEMENT

    /// P1Bubble Element on a 2D triangle
    /*
                           3
                           |\
                           | \
                           |4.\
                           1---2
    */
    class P1Bubble_TriangleElement : public NodalFiniteElement
    {
    public:
        /// Construct the P1Bubble_TriangleElement
        P1Bubble_TriangleElement();
        void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
        void CalcDShape(const IntegrationPoint &ip,
                        DenseMatrix &dshape) const override;
        void CalcHessian(const IntegrationPoint &ip,
                         DenseMatrix &Hessian) const override;
        void ProjectDelta(int vertex, Vector &dofs) const override
        { dofs = 0.0; dofs(vertex) = 1.0; } // { dofs = 1.0; }
    };



    //////////////////////////////////////////////////////////////////////////
    ///                 P1Bubble Element on a 2D triangle                 ///
    //////////////////////////////////////////////////////////////////////////

    // Implementation of the P1Bubble_TriangleElement constructor
    P1Bubble_TriangleElement::P1Bubble_TriangleElement()
        : NodalFiniteElement(2, Geometry::TRIANGLE, 4, 1)
    {
        Nodes.IntPoint(0).x = 0.0;
        Nodes.IntPoint(0).y = 0.0;
        Nodes.IntPoint(1).x = 1.0;
        Nodes.IntPoint(1).y = 0.0;
        Nodes.IntPoint(2).x = 0.0;
        Nodes.IntPoint(2).y = 1.0;
        Nodes.IntPoint(3).x = 0.333333333333333333;
        Nodes.IntPoint(3).y = 0.333333333333333333;
    }

    // Implementation of the P1Bubble_TriangleElement::CalcShape method
    void P1Bubble_TriangleElement::CalcShape(const IntegrationPoint &ip,
                                             Vector &shape) const
    {
        real_t r = ip.x;
        real_t s = ip.y;
        shape(0) = 1.0 - r - s;
        shape(1) = r;
        shape(2) = s;
        // bubble function
        shape(3) = shape(0) * shape(1) * shape(2); // (1 - r - s) * r * s;
    }

    // Implementation of the P1Bubble_TriangleElement::CalcDShape method
    void P1Bubble_TriangleElement::CalcDShape(const IntegrationPoint &ip,
                                              DenseMatrix &dshape) const
    {
        real_t r = ip.x;
        real_t s = ip.y;
        dshape = 0.0;

        dshape(0, 0) = -1.0;                  dshape(0, 1) = -1.0;                // phi1_x, phi1_y
        dshape(1, 0) = 1.0;                   dshape(1, 1) = 0.0;                 // phi2_x, phi2_y
        dshape(2, 0) = 0.0;                   dshape(2, 1) = 1.0;                 // phi3_x, phi3_y
        dshape(3, 0) = (1.0 - 2.0 * r - s)*s; dshape(3, 1) = (1 - r - 2 * s) * r; // phi4_x, phi4_y
    }

    // Implementation of the P1Bubble_TriangleElement::Calch method
    //       The order in 2D is {u_xx, u_xy, u_yy}.
    void P1Bubble_TriangleElement::CalcHessian(const IntegrationPoint &ip,
                                               DenseMatrix &h) const
    {
        real_t r = ip.x;
        real_t s = ip.y;
        h = 0.0;

        h(0, 0) = 0.0;     h(0, 1) = 0.0;                h(0, 2) = 0.0;    // phi1_xx, phi1_xy, phi1_yy
        h(1, 0) = 0.0;     h(1, 1) = 0.0;                h(1, 2) = 0.0;    // phi2_xx, phi2_xy, phi2_yy
        h(2, 0) = 0.0;     h(2, 1) = 0.0;                h(2, 2) = 0.0;    // phi3_xx, phi3_xy, phi3_yy
        h(3, 0) = -2.0*s;  h(3, 1) = 1.0-2.0*r-2.0*s;  h(3, 2) = -2.0*r;   // phi4_xx, phi4_xy, phi4_yy
    }

The FINITE ELEMENT COLLECTION

/// Collection for Linear FE enriched with bubble functions in 2D. (MINI element)
class MINI2D_FECollection : public FiniteElementCollection
{
private:
   const PointFiniteElement PointFE;
   const Linear1DFiniteElement SegmentFE;
   const P1Bubble_TriangleElement TriangleFE;

public:

   MINI2D_FECollection() : FiniteElementCollection(1) {}

   const FiniteElement *
   FiniteElementForGeometry(Geometry::Type GeomType) const override;

   int DofForGeometry(Geometry::Type GeomType) const override;

   const int *DofOrderForOrientation(Geometry::Type GeomType, int Or) const override;

   const char *Name() const override { return "MINI2D"; }

   int GetContType() const override { return CONTINUOUS; }
};




const FiniteElement * MINI2D_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
{
   switch (GeomType)
   {
      case Geometry::POINT:       return &PointFE;
      case Geometry::SEGMENT:     return &SegmentFE;
      case Geometry::TRIANGLE:    return &TriangleFE;
      default:
         if (error_mode == RETURN_NULL) { return nullptr; }
         mfem_error ("MINI2D_FECollection: unknown geometry type.");
   }
   return &TriangleFE; // Make some compilers happy
}

int MINI2D_FECollection::DofForGeometry(Geometry::Type GeomType) const
{
   switch (GeomType)
   {
      case Geometry::POINT:       return 1;
      case Geometry::SEGMENT:     return 0;
      case Geometry::TRIANGLE:    return 1;
      default:
         if (error_mode == RETURN_NULL) { return 0; }
         mfem_error ("MINI2D_FECollection: unknown geometry type.");
   }
   return 1; // Make some compilers happy
}

const int *MINI2D_FECollection::DofOrderForOrientation(Geometry::Type GeomType,
                                                       int Or) const
{
   static int indexes[] = {0};

   return indexes;
}

[tzanio]

ping: @mlstowell, @jandrej