Hello,
Wondering if it is possible to use node values of FEs in varf. Something like the following.
Thank you!
//
mesh Th=square(10,10);
fespace V1h(Th,P1);
V1h u, uh;
varf ma(u,uh)
= int2d(Th) (u[Th[nuTriangle][0]]*uh[u[Th[nuTriangle][0]]]);
Hello,
One possibility is to modify the quadrature formula for triangles as
int2d(Th,qft=qf1pTlump)
With qf1pTlump the integral of f over a triangle is computed by |T|(f(V1)+f(V2)+f(V3))/3
where V_i are the vertices of the triangle.
If you put a cutoff function \phi you can get quite a large number of possibilities writing
int2d(Th,qft=qf1pTlump)(phi*u*uh)
Thank you for the suggestion.
I have thought about it, but it does not solve my problem, because u and uh are on different meshes and the compiler complains.
The meshes are “equivalent” in the sense that one is movemesh of the other. In this case mapu and mapt work, but I prefer to build the matrix in a direct way.
It seems that the modification of the matrix M from varf is the way. But this requires knowing the ordering of dofs, ie to which dof corresponds the entry M(i,j), which in the case of vector FE spaces, such as [P1b,P1b,P1] seems not available in the docs - I have posted a message with this regard at Order of variables in matrix generated with varf with a [P1, P1, P0] fespace - #6 by pi.3.141.
When you write (for example)
fespace X(Th,[P1b,P1b,P1]);
fespace Y(Th,[P1,P1]);
varf aa([x1,x2,x3],[y1,y2]) = int2d(Th)(…);
matrix M=aa(X,Y);
the coefficient M(i,j) of M corresponds to
i (line number) = index of dof of Y
j (column number) = index of dof of X
Then in order to understand the dofs of a product space like X above, you can make a loop over triangles and dofs of X in each triangle, as
for (int k=0;k<Th.nt;++k) {// loop on triangles of the mesh
for (int j=0;j<X.ndofK;++j) {//loop over the dofs in triangle k
cout << "X(k,j)=" << X(k,j) << endl;
}
}
Then X(k,j) is the index of the dof of X (which is required above as column number of M).
The incrementation of j from 0 to X.ndofK-1 corresponds in order to first the dofs in triangle of the first space in product (here P1b), then listing the dofs in triangle of the second space in product (here P1b), then the dofs in triangle of the third space in product (here P1).
The order of dofs in triangle for each space (P1b or P1 or other) is specific to this space.
For P1, the order corresponds to listing the vertices 0,1,2, as numbered by Th[k][j]
For P1b, the order is first listing the 3 vertices 0,1,2 as for P1, then the bubble.
For P2, the order is first the 3 vertices 0,1,2, then the three edge middles, which are in order opposite to the vertices 0,1,2.
An example is here Order of variables in matrix generated with varf with a [P1, P1, P0] fespace - #2 by fb77
I think you have all the information for your purpose.
Thank you, Francois, for the detailed information!
Hello Francois,
I have a similar question related to the composite finite elements.
For example, considering the script in Composite finite element spaces NEW!,
fespace Xh=Uh*Ph;
varf Stokes (<[u1,u2],[p]>, <[v1,v2],[q]>)
= int2d(ThU)((Grad(u1,u2):Grad(v1,v2)))
**+ int2d(ThU)(-div(u1,u2)q -div(v1,v2)p)
+ int2d(ThP)(-1e-10pq)
+ int2d(ThU)([f1,f2]'*[v1,v2])
+ on(1,2,3,4, u1=g1, u2=g2);
matrix M = Stokes(Xh,Xh);
what is the order of dof of Xh?
Thank you!
Further to my question, I am copying a comment from the document
Composite finite element spaces NEW!, which says
“Under the hood, the matrix and right-hand side of a composite problem are assembled in a contiguous way with respect to the different components. This means that the linear system we are solving has the expected structure”.
This helps. However, as in composite FE we are allowed to use quite various spaces, the question is still relevant.
For example, if we have something like
Th=Th1+Th2;
fespace UVh(Th,P1b);
fespace P1h(Th1,P1);
fespace P2h(Th2,P1);
Xh=UVhP1hP2h;
varf stokes(<[u,v],[p1],[p2]>,<[uh,vh],[p1h],[p2h]>)
= int2d(Th)(…) + int2d(Th1)(…) + int2d(Th2)(…) + …;
matrix M=stokes(Xh,Xh);
some more insights are needed to better understand how the dofs are ordered.
I have not tested the order of dofs for composite spaces, but as you say, the document description is quite clear.
I understand that each space in the product has its dofs. Then you take globally the dofs of the first space, then the dofs of the second space, and so on. This is the most simple, and it is the same as when you write a matrix by blocks.
Thank you, Francois!
Following the matter of composite FEs, in this piece of script in the documentation,
1 mesh ThP = square(nn, nn, [2pix,2piy], flags=3);
2 mesh ThU = trunc(ThP,1,split=2);
3 fespace Uh(ThU, [P1,P1]);
4 fespace Ph(ThP, P1);
5 Uh [u,v], [uh,vh];
6 Ph p, ph;
…
7 fespace Xh=UhPh;
8 varf Stokes (<[u,v],[p]>, <[uh,vh],[ph]>)
9 = int2d(ThU)((Grad(u,v):Grad(uh,vh)))
10 + int2d(ThU)(-div(u,v)ph -div(uh,vh)p)
11 + int2d(ThP)(-1e-10pph)
12 - int2d(ThU)([f1,f2]'[uh,vh])
13 + on(1,2,3,4, u=g1, v=g2);
it is not clear the rationale of the choice of the mesh under the integrals.
For example, while the choice of the mesh ThU in line 9 is meaningful, in line 10 it is not clear. One would break the integral of line 10 in two parts, separating (uh,vh) and ph, and taking the respective integrals in ThU and ThP.
This said, taking in line 10, ThU, ThP, or even breaking it in two integrals, one in ThU and the other in ThP, the result is almost the same (velocity is the same while pressure differs by 10^-5 or less).
It would be helpful if you post any comment on this matter.
Thank you!
You are right about the unclear choice of domain of integration. The interest of composite spaces is that you can have together finite element functions defined on different domains. Thus you have a finite element function, for example p, that is defined a priori on ThP, but you want to integrate it on ThU, which is obtained here as ThU=trunc(ThP,1,split=2).
Since the composite spaces were designed exactly in order to do that, you can hope that it does it correctly, which means it evaluates correctly the functions.
However you have to understand that when you write int2d(ThU) (resp int2d(ThP)) the integral is computed as the sum of integrals over the triangles of ThU (resp of ThP). In particular when you integrate on ThP (the coarse mesh) a function which is P1 on ThU (the fine mesh), the result is not necessarily exact since this function is not polynomial on the triangles of ThP. The most accurate is to integrate always over the finer mesh.