Is it possible to access values of a P1 element at the nodes?

General Discussion
1

Hi there!
Is it possible to access values of a P1 function at the nodes and put them in an array?

Thank you.

1 Like
2

Finite element function can be automatically converted in array, you have to write:

fespace Uh(Th, P1);
Uh u = ...;

cout << u[] << endl;
cout << u[][5] << endl;

u[] is the array of u, u[][5] is the 6th element of the array.

2 Likes
3

Thank you for the quick answer.

I have another question. Given a mesh Th and a function u in Vh, how do I get the value of u at the node Th[i][j]?

4

The mesh node is obtained using Th(i) (not Th[i][j], this is triangle i, vertex j, see documentation).

You have access to the value at node i using u[][i], or alternatively u(Th(i).x, Th(i).y) (using interpolation).

1 Like
5

How about in the case of a vector

fespace Vh(Th,[P1,P1]);
Vh [ux, uy];

I gues it should be the same but when I try to output this to file like the following code:

{
            ofstream cupf("file.txt");
            cupf.precision(8);
            cupf.scientific;
            for (int i=0; i < Th.nv; i++)
            {
                cupf << Th(i).x << " " << Th(i).y << " " << ux[][i] << " " << uy[][i] << endl;
            }
            cupf.flush;  //to flush the buffer of file
}

It gives me the exact same values for ux than for uy, which I know is not true because I also plot this vector field and it looks OK. Maybe I’m doing something wrong with the way I access the values since it is a vector? Any ideas are highly appreciated!

6

You have to write:

ux[][2*i] << " " << ux[][2*i+1]

Remark: ux[] and uy[] are the same objects

2 Likes
7

Thanks a lot! It worked. I think this should be in the manual since it is quite confusing that you can define the vectors using component names for the variables but accessing the elements has a different format.

8

Thanks for your feedback.
I will add a clearer description of vectorial finite elements soon in the documentation

1 Like
9

Your welcome and thanks to you.

Another related question, now for P2 elements. In the same script along with the vectors I also have P2 scalar variables:

fespace Vh(Th,[P1,P1]);
fespace Sh2(Th, P2);
Vh [ux, uy];
Sh2 c;

But when I try to write the values to file in the same way and together with the vector values like this, using your suggestion

{
                ofstream cupf("file.txt");
                cupf.precision(16);
                cupf.scientific;
                for (int i=0; i < Th.nv; i++)
                {
                    cupf << Th(i).x << " " << Th(i).y << " " << c[][i] << " " << ux[][2*i] << " " << ux[][2*i+1] << endl;
                }
                cupf.flush;     //to flush the buffer of file
}

The vector component values seem OK now but the scalar value, c variable, does not. The plot directly from FreeFem built in command looks fine but if I plot the scalar c from the data file (I’m working with Matplotlib in Python) it looks like the same values but different distribution, like scrambled…

Figure_2
Figure_2.png929Ă—978 36.6 KB

maybe the order is different because it’s a different type of element? even if the mesh is the same?

10

It’s more complicated with a P2 finite element: your degrees of freedom are not only the mesh nodes.
A simple way to do that is to write:

c(Th(i).x, Th(i).y)

(That use interpolation)

11

It didn’t work. I have the table value in txt file written in format “x \t y \t fxy”, and read the data as per following code:
mesh Th2 = triangulate(tempFileName);
fespace Vh2(Th2,P1);
Vh2 fxy;
{
ifstream temp(tempFileName);
for(int i = 0; i < fxy.n; i++)
{
real x1,y1;
temp >> x1 >> y1 >> fxy[i]; //to read third column only
}
}
Looking at plot of the fxy, the file is read normally and right. Then, I needed function values in vertices of a triangle containing the described coordinate, say, (0, -39):
int N = Th2(0,-39).nuTriangle;
cout << "Number of triangle containing (0,-39) is "<< N << endl;
for (int i = 0; i<3; ++i)
{
int xi = Th2[N][i].x;
int yi = Th2[N][i].y;
cout << “Vertix” << i << "x = " << xi << ", y = " << yi << endl;
cout << "Function value is " << fxy(Th2(N).x, Th2(N).y << endl;
}
I received result
Vertix0x = 1, y = -39
Function value is 0.919753
while the function in stated coordinates (x y fxy) is 1 -39 0.6419752836227417.
I tried also fxy[N] or fxy[i], it gave me N-th and i-th values of function in the table (I mean, values in N-th and i-th row of the table) regardless the vertix coordinate.
Attempt fxy[N][i] failed: Error line number 131, in file c:\projects\FreeFEM\Bubble-1.edp, before token ]. Seems like such syntax is not allowed.

12

But the following appeared as working:
cout << "Function value is " << fxy(Th2[N][i].x,Th2[N][i].y) << endl;