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.
For P1/P0, splitmesh6PowellSabin in 2d, splitmesh12WorseyFarin in 3d. V_h = P1 for velocity, however, P0 is not the true space M_h for pressure. Actually, M_h is the subspace of P0, thus \nable \cdot V_h \subset M_h , which satisfies the divergence-free. So, how to define the true finite element space by using the FreeFem. I read your code given the above. You define the original mesh and then split it. Define the finite element space based on “For P1/P0, splitmesh6PowellSabin in 2d, splitmesh12WorseyFarin in 3d”. Finally, you define the weak formulation for Stokes like using the Taylor-Hood element pair. However, the space P0 is not the real space M_h for pressure. How do you deal with it?
The paper “Low-order divergence-free approximations for the Stokes problem on
Worsey–Farin and Powell–Sabin splits, Comput. Methods Appl. Mech. Engrg. 390 (2022) 114444“ gives an idea for achieving the real M_h. Your above code is consistent with the idea? If not, could you have an idea to achieve this and give me a code example? Thanks for your attention!
Dear Wenyue,
You are right, the case P1/P0 is more involved. Here is the way I do it.
The space M_h for pressure is a subspace of P0, defined locally by some linear relations between the dofs of P0(Th). Th is the mesh obtained by Powell-Sabin split of the mesh Th0.
When solving \int q\mathop{\rm div} u=0, not all q\in P0 are necessary, we need only q\in M_h, since M_h=\mathop{\rm div} V_h. Thus some lines of the matrix of the Stokes problem set on P0 are unnecessary. The principle is to replace these lines by the linear relations that we want to impose on the dofs of the pressure p.
However, implementing this is not so easy because we need to identify the local dofs. They are related to the geometry of the split Th0–>Th. Thus we need to identify these geometrical relations. This is done in the section “pressure space reduction”.
Moreover we have to be careful of having a nonsingular system, and take care of the arbitrary constant in the pressure, and the fact the the test function q=1 is not necessary.
Code: stokes-powell-sabin.edp (12.3 KB)
This code is for the Stokes problem -\mathop{\rm div}Du+\nabla p=f, \qquad\mathop{\rm div}u=0,\qquad with Du\equiv (\nabla u+(\nabla u)^t)/2,
with Dirichlet boundary conditions.
Unfortunately if we change the boundary conditions, for example take Neumann, then the space reduction has to be modified for the boundary edges.
In the code there is also a simplified approach (obtained by taking “macro FEchoice()P1/P0eps” instead of “macro FEchoice()P1/P0”).
The principle is to replace \mathop{\rm div}u=0 by \mathop{\rm div}u+\varepsilon p=0, for some small \varepsilon (taken 1e-8 in the code). Then we can simply get p=-\mathop{\rm div}u/\varepsilon and get the equation on u as -\mathop{\rm div}Du-\frac{1}{\varepsilon}\nabla\mathop{\rm div}u =f.
This method of elimination of pressure is rather efficient, although it does not give exact free divergence. Notice that the Powell-Sabin split is here necessary in order to get the inf-sup condition that ensures the stability of the scheme.
Professor, thank you very much. I ran the code and found that the divergence of the velocity field is very small, about O( 10^(-15) ). So the velocity is exactly div-free. However, how to implement it is also difficult. So, if we use the P2/P1dc with barycenter refinement of the mesh, is it necessary to perform “pressure space reduction“? The code using P2/P1dc to perform numerical tests is similar to the usual Taylor-Hood finite element? Thanks again.
The scheme with pressure space reduction is complicate. This is why I have proposed the simplified approach, which does not need it. Then the div is of order \varepsilon. This parameter has to be small, but not too much otherwise the error in the momentum equation becomes not negligible. I proposed to take a medium value \varepsilon=1e-8=\sqrt{1e-16} in order to balance the two errors (error in the div and error in the momentum). It leads to a divergence error of about 1e-8, which is good enough.
In the case of P2/P1dc with splitmesh3, there is no need of pressure space reduction. The algorithm is simple in this case. stokes-hct.edp (6.5 KB)
the exact case is macro linearsystem()velpressmatrix//
the simplified case (similar to the P1/P0 above) is macro linearsystem()eliminatepressure//
another simplified is macro linearsystem()directsolve// (more intuitive, but most costly than the previous)