About Scott–Vogelius element

lailai · September 10, 2022, 5:09am

Dear all,

Hi, everyone! I have a question regarding the Scott–Vogelius lowest order element for Stokes flows, namely P2 element for velocity and P1dc element (piecewise linear discontinuous Galerkin) for pressure. This element does not satisfy the LBB or inf-sup condition, which can be circumvented using the macroelement technique by splitting a triangular element into three sub triangles. If I understand correctly, the splitting process is available in Freefem as I have seen the examples ‘splitmesh*.edp’ in the /plugin dir. I would like to know whether my above understanding is correct and whether this Scott–Vogelius P2-P1dc element can indeed be adopted in freefem while satisfying the inf-sup condition. Thanks in advance.

best,
lailai

prj · September 10, 2022, 7:02am

Your understanding is correct, it’s perfectly fine to use Scott–Vogelius finite elements and Powel–Sabin mesh refinement. This works in 2D and 3D.

lailai · September 10, 2022, 7:17am

Nice, thank you very much for your prompt response. May I know whether we have a relevant freefem example to learn from? Maybe related to the Powel-Sabin mesh refinement you referred to? Thank you very much.
best,
lailai

prj · September 10, 2022, 7:32am

I’m not sure there are any example, but you can take any Navier–Stokes example, replace a Taylor–Hood pair with a Scott–Vogelius one, and apply Powel–Sabin mesh refinement before discretizing your system.

fb77 · February 2, 2024, 1:51pm

Dear all,
regarding the Powell-Sabin division, I have written some FreeFem lines that work. This can be helpful for some of you. Given a mesh Th0 it produces a mesh Th by subdivision of each initial triangle into 6.
Best,
Francois.

// produce mesh1.msh as a Powell-Sabin division of Th0
// Cf Shangyou Zhang J. Comput. Math. 26 no3, 456-470, 2008.
int[int] dofmesh1vertex0(Th0.nv); //new numbering of old vertices
dofmesh1vertex0=-1;// value -1 means not yet defined
real[int] dofmesh1x(7*Th0.nt),dofmesh1y(7*Th0.nt);// x,y components of new vertices
int[int] dofmesh1label(7*Th0.nt);// labels of new vertices
dofmesh1label=-1;
int[int,int] edgesmesh0(Th0.nt,3);//new numbering of old edges, in terms of old triangles and edges
int[int,int] trianglesmesh1(6*Th0.nt,4);//list of new triangles, by 3 vertices and label
int countdofmesh1=0;
int[int]countcenter(Th0.nt);//new numbering of old triangles
for (int k=0;k<Th0.nt;k++) {//begin loop on triangles
 int v0=Th0[k][0];
 int v1=Th0[k][1];
 int v2=Th0[k][2];
 real v0x=Th0[k][0].x;
 real v0y=Th0[k][0].y;
 real v1x=Th0[k][1].x;
 real v1y=Th0[k][1].y;
 real v2x=Th0[k][2].x;
 real v2y=Th0[k][2].y;
// vertices of Th0
if (dofmesh1vertex0(v0)==-1) {
dofmesh1x(countdofmesh1)=v0x;
dofmesh1y(countdofmesh1)=v0y;
dofmesh1label(countdofmesh1)=Th0[k][0].label;
dofmesh1vertex0(v0)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
}
if (dofmesh1vertex0(v1)==-1) {
dofmesh1x(countdofmesh1)=v1x;
dofmesh1y(countdofmesh1)=v1y;
dofmesh1label(countdofmesh1)=Th0[k][1].label;
dofmesh1vertex0(v1)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
}
if (dofmesh1vertex0(v2)==-1) {
dofmesh1x(countdofmesh1)=v2x;
dofmesh1y(countdofmesh1)=v2y;
dofmesh1label(countdofmesh1)=Th0[k][2].label;
dofmesh1vertex0(v2)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
}
//// center of triangles of Th0
//// gravity center
//dofmesh1x(countdofmesh1)=(v0x+v1x+v2x)/3.;
//dofmesh1y(countdofmesh1)=(v0y+v1y+v2y)/3.;
// center of the inscribed circle
real length0,length1,length2,w0x,w0y,w1x,w1y;
length0=sqrt((v2x-v1x)^2+(v2y-v1y)^2);
length1=sqrt((v0x-v2x)^2+(v0y-v2y)^2);
length2=sqrt((v1x-v0x)^2+(v1y-v0y)^2);
w0x=(v1x-v0x)/length2+(v2x-v0x)/length1;
w0y=(v1y-v0y)/length2+(v2y-v0y)/length1;
w1x=(v2x-v1x)/length0+(v0x-v1x)/length2;
w1y=(v2y-v1y)/length0+(v0y-v1y)/length2;
real lambda0=((v1x-v0x)*w1y-(v1y-v0y)*w1x)/(w0x*w1y-w0y*w1x);
dofmesh1x(countdofmesh1)=v0x+lambda0*w0x;
dofmesh1y(countdofmesh1)=v0y+lambda0*w0y;//end of def center of the inscribed circle
dofmesh1label(countdofmesh1)=Th0[k].label;
countcenter(k)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
// edges of Th0 (cut by line between two centers)
int ee;
if (int(Th0[k].adj((ee=0)))>=k) {
if (int(Th0[k].adj((ee=0)))!=k){
dofmesh1label(countdofmesh1)=0;
}
else {
dofmesh1x(countdofmesh1)=(v1x+v2x)/2.;
dofmesh1y(countdofmesh1)=(v1y+v2y)/2.;
}
edgesmesh0(k,0)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
}
else {
int k0=int(Th0[k].adj((ee=0)));
edgesmesh0(k,0)=edgesmesh0(k0,ee);
real centerkx0=dofmesh1x(countcenter(k));
real centerky0=dofmesh1y(countcenter(k));
real centeradjx0=dofmesh1x(countcenter(k0));
real centeradjy0=dofmesh1y(countcenter(k0));
real scalinters0=((v1x-centerkx0)*(v2y-v1y)-(v1y-centerky0)*(v2x-v1x))/((centeradjx0-centerkx0)*(v2y-v1y)-(centeradjy0-centerky0)*(v2x-v1x));
dofmesh1x(edgesmesh0(k0,ee))=centerkx0+scalinters0*(centeradjx0-centerkx0);
dofmesh1y(edgesmesh0(k0,ee))=centerky0+scalinters0*(centeradjy0-centerky0);
}
if (int(Th0[k].adj((ee=1)))>=k) {
if (int(Th0[k].adj((ee=1)))!=k){
dofmesh1label(countdofmesh1)=0;
}
else {
dofmesh1x(countdofmesh1)=(v2x+v0x)/2.;
dofmesh1y(countdofmesh1)=(v2y+v0y)/2.;
}
edgesmesh0(k,1)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
}
else {
int k1=int(Th0[k].adj((ee=1)));
edgesmesh0(k,1)=edgesmesh0(k1,ee);
real centerkx1=dofmesh1x(countcenter(k));
real centerky1=dofmesh1y(countcenter(k));
real centeradjx1=dofmesh1x(countcenter(k1));
real centeradjy1=dofmesh1y(countcenter(k1));
real scalinters1=((v2x-centerkx1)*(v0y-v2y)-(v2y-centerky1)*(v0x-v2x))/((centeradjx1-centerkx1)*(v0y-v2y)-(centeradjy1-centerky1)*(v0x-v2x));
dofmesh1x(edgesmesh0(k1,ee))=centerkx1+scalinters1*(centeradjx1-centerkx1);
dofmesh1y(edgesmesh0(k1,ee))=centerky1+scalinters1*(centeradjy1-centerky1);
}
if (int(Th0[k].adj((ee=2)))>=k) {
if (int(Th0[k].adj((ee=2)))!=k){
dofmesh1label(countdofmesh1)=0;
}
else {
dofmesh1x(countdofmesh1)=(v0x+v1x)/2.;
dofmesh1y(countdofmesh1)=(v0y+v1y)/2.;
}
edgesmesh0(k,2)=countdofmesh1;
countdofmesh1=countdofmesh1+1;
}
else {
int k2=int(Th0[k].adj((ee=2)));
edgesmesh0(k,2)=edgesmesh0(k2,ee);
real centerkx2=dofmesh1x(countcenter(k));
real centerky2=dofmesh1y(countcenter(k));
real centeradjx2=dofmesh1x(countcenter(k2));
real centeradjy2=dofmesh1y(countcenter(k2));
real scalinters2=((v0x-centerkx2)*(v1y-v0y)-(v0y-centerky2)*(v1x-v0x))/((centeradjx2-centerkx2)*(v1y-v0y)-(centeradjy2-centerky2)*(v1x-v0x));
dofmesh1x(edgesmesh0(k2,ee))=centerkx2+scalinters2*(centeradjx2-centerkx2);
dofmesh1y(edgesmesh0(k2,ee))=centerky2+scalinters2*(centeradjy2-centerky2);
}
}//end of loop on triangles
//cout << "countdofmesh1= " << countdofmesh1 << endl;
dofmesh1x.resize(countdofmesh1);
dofmesh1y.resize(countdofmesh1);
dofmesh1label.resize(countdofmesh1);

// def of 6 subtriangles
for (int k=0;k<Th0.nt;k++) {
trianglesmesh1(6*k,0)=dofmesh1vertex0(Th0[k][0]);
trianglesmesh1(6*k,1)=edgesmesh0(k,2);
trianglesmesh1(6*k,2)=countcenter(k);
trianglesmesh1(6*k,3)=Th0[k].label;
trianglesmesh1(6*k+1,0)=dofmesh1vertex0(Th0[k][1]);
trianglesmesh1(6*k+1,1)=countcenter(k);
trianglesmesh1(6*k+1,2)=edgesmesh0(k,2);
trianglesmesh1(6*k+1,3)=Th0[k].label;
trianglesmesh1(6*k+2,0)=dofmesh1vertex0(Th0[k][1]);
trianglesmesh1(6*k+2,1)=edgesmesh0(k,0);
trianglesmesh1(6*k+2,2)=countcenter(k);
trianglesmesh1(6*k+2,3)=Th0[k].label;
trianglesmesh1(6*k+3,0)=dofmesh1vertex0(Th0[k][2]);
trianglesmesh1(6*k+3,1)=countcenter(k);
trianglesmesh1(6*k+3,2)=edgesmesh0(k,0);
trianglesmesh1(6*k+3,3)=Th0[k].label;
trianglesmesh1(6*k+4,0)=dofmesh1vertex0(Th0[k][2]);
trianglesmesh1(6*k+4,1)=edgesmesh0(k,1);
trianglesmesh1(6*k+4,2)=countcenter(k);
trianglesmesh1(6*k+4,3)=Th0[k].label;
trianglesmesh1(6*k+5,0)=dofmesh1vertex0(Th0[k][0]);
trianglesmesh1(6*k+5,1)=countcenter(k);
trianglesmesh1(6*k+5,2)=edgesmesh0(k,1);
trianglesmesh1(6*k+5,3)=Th0[k].label;
}
int[int,int]vertexbdryel(2*Th0.nbe,3);//list of new boundary elements, by 2 vertices and label
// loop on boundary
for (int l=0;l<Th0.nbe;l++) {
 vertexbdryel(2*l,0)=dofmesh1vertex0(Th0.be(l)[0]);
 int triel=int(Th0.be(l).Element);
 int edgetriel=Th0.be(l).whoinElement;
 vertexbdryel(2*l,1)=edgesmesh0(triel,edgetriel);
 vertexbdryel(2*l,2)=Th0.be(l).label;
 vertexbdryel(2*l+1,0)=edgesmesh0(triel,edgetriel);
 vertexbdryel(2*l+1,1)=dofmesh1vertex0(Th0.be(l)[1]);
 vertexbdryel(2*l+1,2)=Th0.be(l).label;
 dofmesh1label(edgesmesh0(triel,edgetriel))=Th0.be(l).label;
}

{ofstream filemesh1("mesh1.msh");
 filemesh1 << countdofmesh1 << " " << 6*Th0.nt << " " << 2*Th0.nbe << endl;
};
{ofstream filemesh1("mesh1.msh",append);
int nold=filemesh1.precision(30);// it is important to have really straight intersections
for (int i=0;i<countdofmesh1;i++) {
 filemesh1 << dofmesh1x(i) << " " << dofmesh1y(i) << " " << dofmesh1label(i) << endl;
}
for (int k=0;k<6*Th0.nt;k++) {
 filemesh1 << trianglesmesh1(k,0)+1 << " " << trianglesmesh1(k,1)+1 << " " << trianglesmesh1(k,2)+1 << " " << trianglesmesh1(k,3) << endl;
}
for (int l=0;l<2*Th0.nbe;l++) {
 filemesh1 << vertexbdryel(l,0)+1 << " " << vertexbdryel(l,1)+1 << " " << vertexbdryel(l,2) << endl;
}
};
//end of Powell-Sabin division

mesh Th=readmesh("mesh1.msh");

prj · February 2, 2024, 2:05pm

Why not use the splitmesh6 plugin? It’s exactly what it does, but it’s order of magnitude faster than your code. Cf. examples/plugin/splitmesh6.edp.

fb77 · February 2, 2024, 2:08pm

Thank you. I do not know about the content of splitmesh6. Do you have some description of what it does exactly ? Or is there a source of it (in FreeFem ?)

prj · February 2, 2024, 2:13pm

All splitmesh* codes are for computing Powell–Sabin mesh refinements. See:

fb77 · February 2, 2024, 2:22pm

This gives examples of use, but no description of what it does nor any source code.

prj · February 2, 2024, 2:23pm

I gave you the description, the code is in the plugin folder, as all other plugins…

fb77 · February 2, 2024, 2:34pm

Sorry, I’m not very comfortable with the Github. I do not see the source of the splitmesh6 procedure, but only the code of the example of use.

prj · February 2, 2024, 2:36pm

github.com

FreeFem/FreeFem-sources/blob/220957ed86797a6c1ed771f1fe1c27078e94f28a/plugin/seq/splitmesh6.cpp

/****************************************************************************/
/* This file is part of FreeFEM.                                            */
/*                                                                          */
/* FreeFEM is free software: you can redistribute it and/or modify          */
/* it under the terms of the GNU Lesser General Public License as           */
/* published by the Free Software Foundation, either version 3 of           */
/* the License, or (at your option) any later version.                      */
/*                                                                          */
/* FreeFEM is distributed in the hope that it will be useful,               */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of           */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */
/* GNU Lesser General Public License for more details.                      */
/*                                                                          */
/* You should have received a copy of the GNU Lesser General Public License */
/* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>.          */
/****************************************************************************/
// SUMMARY : ...
// LICENSE : LGPLv3
// ORG     : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
// AUTHORS : ...

This file has been truncated. show original

fb77 · February 2, 2024, 2:49pm

Thanks a lot!
Reading the file I see that this is not the true Powell-Sabin refinement, but the simplified one:
– Center of triangle is taken to be the gravity center
– Points on edges are taken as middle of edges
The one that I have implemented is as follows:
– Center of triangle is taken the center of the inscribed circle
– Points on edges are obtained by intersecting the edge with the line joining the centers of the two adjacent triangles

You posted the result of a simple mesh Th0 made from a Cartesian mesh. On this mesh both procedures give the same refinement.

On general meshes the two procedures do not give the same, and this makes a lot of difference:
– the simplified one fails to satisfy the inf-sup condition on P1/P0 elements for velocity/pressure
on the Stokes problem
– the true one satisfies the inf-sup condition

Francois.

prj · February 2, 2024, 3:01pm

I see, what about Scott–Vogelius pairs?

fb77 · February 2, 2024, 3:21pm

Here is an example: original mesh, simplifiedPS, truePS

On the simplified PS mesh you see some vertices for which the associated edges around it form almost two straight lines, but not exactly (for example on the lower right). This causes the infsup constant to be small, and the failure of Stokes computations with P1/P0.

On the true PS you have really straight lines, then the infsup constant is not small, and the computations succeed for P1/P0.

fb77 · February 2, 2024, 9:07pm

To compare the two meshes on the Scott-Vogelius pair P1/P0, consider the Stokes problem on an arbitrary mesh

// Stokes test

border gama11(t= 0.,  1.){ x= t   ; y = 0.;label=1;};
border gama22(t= 0.,  1.){ x= 1.  ; y = t  ;label=2;};
border gama33(t= 1. , 0.){ x= t   ; y = 1. ;label=3;};
border gama44(t= 1. , 0.){ x= 0.; y = t ;label=4;};

int n=30;

mesh Th0=buildmesh(gama11(n)+gama22(n)+gama33(n)+gama44(n));
plot(Th0,fill=0, ps = "maillage0.eps", wait=1);


// here you have to define Th in terms of Th0
// either by       mesh Th=splitmesh6(Th0);// simplified PS
// or by   the code I posted earlier   // true PS


fespace Uh(Th,P1); Uh u,v,uu,vv;
fespace Ph(Th,P0); Ph p,pp;

solve stokes([u,v,p],[uu,vv,pp]) =
int2d(Th)(dx(u)* dx(uu)+dy(u)* dy(uu) + dx(v)* dx(vv)+ dy(v)* dy(vv)
- p*( dx(uu)+dy(vv)) - pp*(dx(u)+dy(v))
- 1e-10*p*pp)
+ on(1,2,4,u=0,v=0) + on(3,u=1,v=0);


plot(u, ps = "u_solution.eps",wait =1,fill=1,value =1);
plot(v, ps = "v_solution.eps",wait =1,fill=1,value =1);
plot([u,v],wait=1);
plot(p,wait=1,fill=1,value =1);

For the simplified you get a wrong solution, and for the true PS you get the right one

frederichecht · February 4, 2024, 10:05am

Remark, it is very simple to add a flag in splitmesh6 to build the Powell-Sabin division.
In can do is rapidly.

fb77 · February 4, 2024, 12:41pm

This would be great, thanks a lot Frédéric.
François.

frederichecht · February 6, 2024, 4:03pm

I do a first version, but I am not sure on the way to split boundary edge

The one that I have implemented is as follows:
– Center of triangle is taken the center of the inscribed circle
– Points on edges are obtained by intersecting the edge with the line joining the centers of the two

splitmesh6.cpp.zip (3.3 KB)
splitmesh6.edp (768 Bytes)

fb77 · February 6, 2024, 4:51pm

For the boundary edges the splitting is not important, and taking the middle of the boundary edge is the most simple.

Topic		Replies	Views
A higher BDM element in Freefem General Discussion	4	245	June 10, 2023
Implementing new finite element in FreeFEM General Discussion	11	429	March 25, 2022
Poiseuille profil, Stokes , Uzawa General Discussion	4	804	May 25, 2024
Inviscid Incompressible Flow in FreeFEM General Discussion	72	356	June 20, 2024
How to solve the following PDEs in Freefem General Discussion	1	94	June 20, 2024

About Scott–Vogelius element

Related topics