FreeFem / FreeFem-sources

FreeFEM source code
https://freefem.org/
Other
796 stars 193 forks source link

Unexpected behavior with global x, y, z variables for finite element function value #253

Closed sgarnotel closed 1 year ago

sgarnotel commented 1 year ago

I noticed a strange behavior using global x, y and z global variable to get finite element function value at some points.

mesh Th = square(10, 10);
fespace Uh(Th, P1);
Uh u = x*y;

x = 0.5;
y = 0.5;

cout << u << endl;
// Good result
// 0.25

cout << u(x, y) << endl;
// Good result
// 0.25

{
    load "iovtk"
    savevtk("test.vtu", Th, u); // This line break the code
}

x = 0.5;
y = 0.5;
cout << u << endl;
// Wrong result
// 0.033333

cout << u(x, y) << endl;
// Good result
// 0.25

Before I use savevtk, the result is correct using the global variables x and y, equally using u(x, y). But after use savevtk the result is wrong using global variables x and y, but correct using u(x, y)

frederichecht commented 1 year ago

Salut, j'ai trouvé

voila la correction

--- a/plugin/seq/iovtk.cpp +++ b/plugin/seq/iovtk.cpp @@ -1847,7 +1847,7 @@ class VTK_WriteMesh_Op : public E_F0mps {

 void writesolutionP0_float_binary(FILE *fp, const Mesh &Th, Stack stack, bool surface,
                                   bool bigEndian) const {

je vais voir dans les autres fonctions !!!!

Le 30 janv. 2023 à 09:35, Simon Garnotel @.***> a écrit :

I noticed a strange behavior using global x, y and z global variable to get finite element function value at some points.

mesh Th = square(10, 10); fespace Uh(Th, P1); Uh u = x*y;

x = 0.5; y = 0.5;

cout << u << endl; // Good result // 0.25

cout << u(x, y) << endl; // Good result // 0.25

{ load "iovtk" savevtk("test.vtu", Th, u); // This line break the code }

x = 0.5; y = 0.5; cout << u << endl; // Wrong result // 0.033333

cout << u(x, y) << endl; // Good result // 0.25 Before I use savevtk, the result is correct using the global variables x and y, equally using u(x, y). But after use savevtk the result is wrong using global variables x and y, but correct using u(x, y)

— Reply to this email directly, view it on GitHub https://github.com/FreeFem/FreeFem-sources/issues/253, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABLPNDGHWSB2LJ6FKNF3XA3WU54OZANCNFSM6AAAAAAUK3IFDA. You are receiving this because you were assigned.

sgarnotel commented 1 year ago

Just copy/paste to highlight the diff

--- a/plugin/seq/iovtk.cpp
+++ b/plugin/seq/iovtk.cpp
@@ -1847,7 +1847,7 @@ class VTK_WriteMesh_Op : public E_F0mps {

     void writesolutionP0_float_binary(FILE *fp, const Mesh &Th, Stack stack, bool surface,
                                       bool bigEndian) const {
-      MeshPoint *mp3(MeshPointStack(stack));
+      MeshPoint *mp3(MeshPointStack(stack)),mps=*mp3;
       R2 Cdg_hat = R2(1. / 3., 1. / 3.);

       for (int it = 0; it < Th.nt; it++) {
@@ -1886,7 +1886,7 @@ class VTK_WriteMesh_Op : public E_F0mps {
       }

       fprintf(fp, "\n");
-    }
+      *mp3 = mps;}

     void writesolutionP0_float(FILE *fp, const Mesh &Th, Stack stack, bool surface) const {
       MeshPoint *mp3(MeshPointStack(stack));
@@ -2770,6 +2770,8 @@ void VTK_WRITE_MESH(const string &filename, FILE *fp, const Mesh &Th, bool