Running your new plugin I observe that the Powell-Sabin mesh is correct, but there is some difficulty:

When solving the Stokes linear system it takes 6 times more in cpu to get the solution. This indicates a stiff problem.
I think that making the computations of the intersection between the edge and the line joining the two triangle centers in C++ is not accurate enough compared to the same computation in FreeFem++ (that I did).
The consequence is that the points are not accurately enough aligned.
I had such problems before when saving my mesh into mesh1.msh:
I had to indicate a large number of digits when writing into the file.
Is it possible to use higher accuracy in your splitmesh6.cpp script ?
Francois.

Dear teacher, the command “splitmesh6” does produce the powell-sabin mesh (in most situations), but it only suits two dimensions. Is there a similar command for producing the powell-sabin mesh in a three-dimensional area?

In 3d this is called Worsey-Farin, see
Low-order divergence-free approximations for the Stokes problem on
Worsey–Farin and Powell–Sabin splits, Comput. Methods Appl. Mech. Engrg. 390 (2022) 114444

Dear teacher, I am very grateful for your reply. And I get your idea. But when I use the Scott-V element (e.g., P2-P1disc), I need some restrictions on mesh for the inf_sup condition. In 2D, we need the barycenter refinement and the powell-sabin refinement for 3D. So, using W-F mesh (splitmesh12) can also satisfy the LBB condition for the SV element ? Thank you again honestly.

For short, the adapted split to satisfy the inf-sup condition is

– For P1/P0,
splitmesh6PowellSabin in 2d
splitmesh12WorseyFarin in 3d
– For P2/P1dc
splitmesh3 in 2d
splitmesh4 in 3d

Indeed the inf-sup condition is satisfied on V_h\times M_h with V_h either P1 or P2 according to the above, and M_h a certain subspace of P0 or P1dc, containing the image of V_h by the div operator. The precise definition of M_h is not important (it is not the full space), what is important is that it contains \mathop{\rm div}V_h.