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;}
|