# About Scott–Vogelius element

**URL:** <https://community.freefem.org/t/about-scott-vogelius-element/2000>\
**Category:** General Discussion\
**Created:** [September 10, 2022, 5:09am UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000 "2022-09-10T05:09:53Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![lailai](https://avatars.discourse-cdn.com/v4/letter/l/e5b9ba/32.png) [@lailai](https://community.freefem.org/u/lailai)\
**Post date:** [September 10, 2022, 5:09am UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/1 "2022-09-10T05:09:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 10, 2022, 7:02am UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/2 "2022-09-10T07:02:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![lailai](https://avatars.discourse-cdn.com/v4/letter/l/e5b9ba/32.png) [@lailai](https://community.freefem.org/u/lailai)\
**Post date:** [September 10, 2022, 7:17am UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/3 "2022-09-10T07:17:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 10, 2022, 7:32am UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/4 "2022-09-10T07:32:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 1:51pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/5 "2024-02-02T13:51:36Z")

</div>

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.

```auto
// 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");

```

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 2, 2024, 2:05pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/6 "2024-02-02T14:05:27Z")

</div>

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`.

 ![Screenshot 2024-02-02 at 3.06.45 PM](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/a/a4f94254e9d85cdf5d0ced528d32019013bceac8.jpeg) ![Screenshot 2024-02-02 at 3.06.47 PM](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/64085818fe05426e1590eeef0c09a9b463579c2d.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 2:08pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/7 "2024-02-02T14:08:30Z")

</div>

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 ?)

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 2, 2024, 2:13pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/8 "2024-02-02T14:13:51Z")

</div>

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

- [FreeFem-sources/examples/plugin/splitmesh3.edp at 220957ed86797a6c1ed771f1fe1c27078e94f28a · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/220957ed86797a6c1ed771f1fe1c27078e94f28a/examples/plugin/splitmesh3.edp)
- [FreeFem-sources/examples/plugin/splitmesh6.edp at 220957ed86797a6c1ed771f1fe1c27078e94f28a · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/220957ed86797a6c1ed771f1fe1c27078e94f28a/examples/plugin/splitmesh6.edp)
- [FreeFem-sources/examples/plugin/splitmesh4.edp at 220957ed86797a6c1ed771f1fe1c27078e94f28a · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/220957ed86797a6c1ed771f1fe1c27078e94f28a/examples/plugin/splitmesh4.edp)
- [FreeFem-sources/examples/plugin/splitmesh12.edp at 220957ed86797a6c1ed771f1fe1c27078e94f28a · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/220957ed86797a6c1ed771f1fe1c27078e94f28a/examples/plugin/splitmesh12.edp).

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 2:22pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/9 "2024-02-02T14:22:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 2, 2024, 2:23pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/10 "2024-02-02T14:23:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 2:34pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/11 "2024-02-02T14:34:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 2, 2024, 2:36pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/12 "2024-02-02T14:36:12Z")

</div>

> <https://github.com/FreeFem/FreeFem-sources/blob/220957ed86797a6c1ed771f1fe1c27078e94f28a/plugin/seq/splitmesh6.cpp>

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 2:49pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/13 "2024-02-02T14:49:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 2, 2024, 3:01pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/14 "2024-02-02T15:01:47Z")

</div>

I see, what about Scott–Vogelius pairs?

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 3:21pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/15 "2024-02-02T15:21:15Z")

</div>

Here is an example: original mesh, simplifiedPS, truePS

 ![maillage_orig](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/6053bb27ef9138aac87c954c99958b1bf2987330.jpeg)  
 ![maillage_PSsimpl](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/0/0b63f6b66ead6b8367ef3adfb51a0aec35984f49.jpeg)  
 ![maillage_PS](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/5/53b262c967dc076eb33a6b06631b8ce01c24654c.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 2, 2024, 9:07pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/16 "2024-02-02T21:07:01Z")

</div>

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

```auto
// 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

 ![u_PS-simpl](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/3a68fb224be24f30e213344bcc112e92dbaeea1c.jpeg)

 ![u_PS](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/d/dae19e6fdbae86d032cdf4c223c4cbb848673d7f.jpeg)

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [February 4, 2024, 10:05am UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/17 "2024-02-04T10:05:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 4, 2024, 12:41pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/18 "2024-02-04T12:41:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [February 6, 2024, 4:03pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/19 "2024-02-06T16:03:19Z")

</div>

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](https://community.freefem.org/uploads/short-url/tybbhrWV63jJt4lOpNIZ0omNuCo.zip) (3.3 KB)  
[splitmesh6.edp](https://community.freefem.org/uploads/short-url/1q0gU8UnRDyUgpktJfFAp9ecqWP.edp) (768 Bytes)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [February 6, 2024, 4:51pm UTC](https://community.freefem.org/t/about-scott-vogelius-element/2000/20 "2024-02-06T16:51:33Z")

</div>

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

[Next page](https://community.freefem.org/t/about-scott-vogelius-element/2000.md?page=2)