# mesh of a right paving stone and a cube

**URL:** <https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394>\
**Category:** General Discussion\
**Created:** [April 13, 2023, 12:04pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394 "2023-04-13T12:04:44Z")\
**Posts on this page:** 20\
**Page:** 2

<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:** [May 3, 2023, 9:18am UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/21 "2023-05-03T09:18:01Z")

</div>

One idea is to introduce a auxiliary variable c = ||v||  
when you can see the problem as minimisation under constraint and use IPOPT to solve the non  
linear problem for fully implicite schema in time

but otherwise you can explicite c a do a semi implicite schema in time.

1. start you 1d or 2d cas to see if the scheme converge (no idea)
2. try do solve a problem with no dependance in time,

---

<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:** [May 3, 2023, 10:02am UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/22 "2023-05-03T10:02:38Z")

</div>

Juste a small example

[aweni.edp](https://community.freefem.org/uploads/short-url/eQ4IfldsYDTFRNmmaeJ4PfQuzXm.edp) (582 Bytes)

with a fixed point iteration,  
remark, the you can do Newton iteration, in 10 mn.

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 3, 2023, 11:52am UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/23 "2023-05-03T11:52:28Z")

</div>

Thank you very much Professor  
I would like to know what this part of the code is for so that I can do it for the 2d  
case before moving on to the time dependent case.

up-= u;  
real err= up.linfty;  
cout \<\< iter \<\< " c "\<\< c\<\< " err = "\<\< err\<\< endl;  
if(err \< 1e-6) break;

---

<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:** [May 3, 2023, 12:38pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/24 "2023-05-03T12:38:08Z")

</div>

```auto
// here up = previous u 
up-= u; // in up = up -u 
real err= up.linfty; // nom of the diff up-u ..
cout << iter << " c "<< c<< " err = "<< err<< endl;
if(err < 1e-6) break;

```

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 3, 2023, 12:44pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/25 "2023-05-03T12:44:14Z")

</div>

Sorry for the 3d i want to write.  
I would like to understand better so that I don’t have to deal with this kind of problem anymore.

up-= u;  
real err= up.linfty;  
cout \<\< iter \<\< " c "\<\< c\<\< " err = "\<\< err\<\< endl;  
if(err \< 1e-6) break;

---

<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:** [May 3, 2023, 12:58pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/26 "2023-05-03T12:58:00Z")

</div>

the 3d version  
[aweni.edp](https://community.freefem.org/uploads/short-url/b18Fmu3YAkKZGQKSBZuPQMiulai.edp) (605 Bytes)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 3, 2023, 12:58pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/27 "2023-05-03T12:58:59Z")

</div>

AFAIK, linifinty is just a fancy way of saying maximum absolute error  
but you can presumably pick a criteria you like.

The static case is intersting ( or would be interesting )f the signs  
change ) as the period is related to the amplitude. Wavelengths  
tend to be fixed by geometry and boundary conditions too.

Since you speciffied initial conditions over omega you already have an  
expression for the time derivative which is probably not entirely zero.

You could write the thing as v=v\_0 + dt\*( \laplacian v + cv)  
with c derived from v\_0 and hopefully in the limit of small dt results make sense  
of as I have done do a 2d+1 case where the “z” dimension in the mesh is actually time.

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 3, 2023, 2:16pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/28 "2023-05-03T14:16:31Z")

</div>

this is what I wrote and what I obtain justifies well what I obtained in the theory :

mesh Th=square(10,10);  
// solve partial\_t u -Delta u + c u = f , l = || u||\_2  
// u = g on gamma  
// arg min int2d(Th) 1/2 grad(u) grad(u) + 1/2 int2d(Th) u^2  
macro grad(u) [dx(u),dy(u)] ) //  
func g = 0;  
func f= 0.;  
// solution u = 1  
int kk=0;  
real c = 0, dt= 0.1,Tfinal=10;  
real[int] instT(floor(Tfinal/dt)+1);// vecteur qui contiendra les instants de temps tn  
real[int] NORML2(floor(Tfinal/dt)+1);// vecteur qui contiendra les différents valeurs de E(tn)  
fespace Vh(Th,P1);

Vh u=x^2+y^2+1 ,up,v ;  
for (real t=0;t\<=Tfinal ;t+=dt)// boucle du temps t  
{

for (int iter= 0; iter \< 100; ++iter)  
{  
up= u;  
c = sqrt(int2d(Th)(u\*u));  
solve Pb(u,v) = int2d(Th) ((1/dt)_u_v)

- int2d(Th) ( [dx(u),dy(u)]'\* [dx(v),dy(v)] - c_u_v)

- int2d(Th) ((1/dt)_up_v) - int2d(Th)(f_v) + on(1,2,3,4, u=g);  
up[]-= u[];  
real err= up[].linfty;  
cout \<\< iter \<\< " c "\<\< c\<\< " err = "\<\< err\<\< endl;  
if(err \< 1e-6) break;  
plot(u, wait=1,dim=3);  
}  
//Calcul of E(t)= |u|\_L^2(Omega)+|nabla u|\_L^2(Omega)  
instT[kk]= t;  
NORML2[kk] = int2d(Th)( u_u)

- int2d(Th)([dx(u),dy(u)]'\* [dx(u),dy(u)]);  
cout \<\<“t” \<\< instT(kk)\<\< endl;

cout \<\<"E(t) " \<\< NORML2(kk) \<\< endl;  
//save of E(t) and t

if ( !(kk % 10)) // kk%10= le reste de la division kk/10  
{  
{  
ofstream fout(“energie.txt”);  
for(int i=0;i \<= kk ; i++)  
fout \<\< instT(i)\<\<" " \<\< NORML2(i)\<\<endl ;  
}  
}  
kk+=1 ;

}

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 3, 2023, 2:24pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/29 "2023-05-03T14:24:06Z")

</div>

You could write the thing as v=v\_0 + dt\*( \laplacian v + cv)  
with c derived from v\_0 and hopefully in the limit of small dt results make sense  
of as I have done do a 2d+1 case where the “z” dimension in the mesh is actually time.

I don’t understand it is the L^2 norm of v .  
I don’t see how it is the derivative of v0 and the derivative is with respect to which variable

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 3, 2023, 4:59pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/30 "2023-05-03T16:59:11Z")

</div>

The L2 norm is in “c” which I thought was used above. where you had a time  
derivaive. Let me look at the code you just posted when I get a chance.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 3, 2023, 6:30pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/31 "2023-05-03T18:30:53Z")

</div>

Can you post your code as an attached file? Copy/paste makes a lot ofthings to clean up.  
Thanks.  
What does it do in the limit of dt → zero?

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 3, 2023, 11:06pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/32 "2023-05-03T23:06:07Z")

</div>

I will send you the file.  
My objective was to see that the energy |v\|\_L^2 decreases towards 0 when the time t tends towards infinity.

(Attachment ss is missing)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 3, 2023, 11:57pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/33 "2023-05-03T23:57:31Z")

</div>

It says attachment is missing or something similar.

Typically you can test your discretization by letting the steps get smaller until  
the results don’t change. Some of these things are unstable with steps too  
big also.

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 4, 2023, 7:36am UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/34 "2023-05-04T07:36:24Z")

</div>

ok sorry. I will resend it

[ss.edp](https://community.freefem.org/uploads/short-url/fDDmcMJYDCav6xmWswUymNyZjqK.edp) (1.57 KB)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 4, 2023, 10:26am UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/35 "2023-05-04T10:26:48Z")

</div>

I had to edit a lot of your upload so I did not do a diff but there are a number of  
issues. I used one of the idp files to move your parameters to the command line,  
see the uploaded fie.  
First, you appear to be using “up” for convergence testing and time stepping.  
Then, dt is in the wrong place. and you can see the huge difference  
by changing the time step although you expect a lot of this when  
the time step is too big anywya I thought your initial conditions were given and you  
appear to iterate those away. I’ll see if I can change it a bit if I get time later.

I always end up with simple algebra mistakes.  
The integration by parts for the weak form is almost magic- it seems  
to remove a lot of terms when more things vary spatially- but it is also  
another source of sign mistakes in algebra 🙂

```auto
 FreeFem++ forum7.edp | grep "t=" | grep " E(t) " | more
   43 : cout <<"t=" << instT(kk)<<" E(t) " << NORML2(kk) << endl;
t=0 E(t) 0.034371
t=0.001 E(t) 0.00063618
t=0.002 E(t) 1.16154e-05
t=0.003 E(t) 2.11678e-07
t=0.004 E(t) 3.85664e-09
t=0.005 E(t) 5.8702e-10
t=0.006 E(t) 5.63973e-10
t=0.007 E(t) 5.4183e-10
t=0.008 E(t) 5.20557e-10
t=0.009 E(t) 5.00119e-10
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ forum7.edp -dt .0005 | grep "t=" | grep " E(t) " | more
   43 : cout <<"t=" << instT(kk)<<" E(t) " << NORML2(kk) << endl;
t=0 E(t) 0.244795
t=0.0005 E(t) 0.0336824
t=0.001 E(t) 0.00455287
t=0.0015 E(t) 0.000610981
t=0.002 E(t) 8.1772e-05
t=0.0025 E(t) 1.09334e-05
t=0.003 E(t) 1.46133e-06
t=0.0035 E(t) 1.95292e-07
t=0.004 E(t) 2.60976e-08
t=0.0045 E(t) 3.48745e-09
t=0.005 E(t) 2.33179e-09
t=0.0055 E(t) 2.28533e-09
t=0.006 E(t) 2.23979e-09

```

[forum7.edp](https://community.freefem.org/uploads/short-url/nVkfDqfkEySGFxD2NwyOZqjVGr7.edp) (1.7 KB)

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 4, 2023, 11:16am UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/36 "2023-05-04T11:16:44Z")

</div>

What do you think of this ( note that this is probably one or more sign mistakes  
but if you have a theory you can compare it to that lol). Note  
that the result is not too sensitive to dt although ti may still be wrong.  
( note I’m just doing this for practice and any feedback would help ).

```auto
FreeFem++ forum7x.edp -dt .0005 | grep "t=" | grep " E(t) " | tail
t=0.095 E(t) 0.0322261
t=0.0955 E(t) 0.0315783
t=0.096 E(t) 0.0309436
t=0.0965 E(t) 0.0303217
t=0.097 E(t) 0.0297124
t=0.0975 E(t) 0.0291153
t=0.098 E(t) 0.0285303
t=0.0985 E(t) 0.0279571
t=0.099 E(t) 0.0273955
t=0.0995 E(t) 0.0268452
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ forum7x.edp -dt .0002 | grep "t=" | grep " E(t) " | tail
t=0.098 E(t) 0.0285142
t=0.0982 E(t) 0.028283
t=0.0984 E(t) 0.0280536
t=0.0986 E(t) 0.0278261
t=0.0988 E(t) 0.0276005
t=0.099 E(t) 0.0273767
t=0.0992 E(t) 0.0271547
t=0.0994 E(t) 0.0269345
t=0.0996 E(t) 0.0267161
t=0.0998 E(t) 0.0264995
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ forum7x.edp -dt .001 | grep "t=" | grep " E(t) " | tail
t=0.09 E(t) 0.0394554
t=0.091 E(t) 0.0378915
t=0.092 E(t) 0.0363899
t=0.093 E(t) 0.034948
t=0.094 E(t) 0.0335636
t=0.095 E(t) 0.0322342
t=0.096 E(t) 0.0309577
t=0.097 E(t) 0.029732
t=0.098 E(t) 0.0285551
t=0.099 E(t) 0.0274249

```

[forum7x.edp](https://community.freefem.org/uploads/short-url/3HZfbHAxfyyG8hA0CPrUwywC7La.edp) (2.0 KB)

```auto
 diff forum7*
11c11
< macro grad(u) [dx(u),dy(u)] ) //
---
> macro grad(u) [dx(u),dy(u)] //
17c17
< real Tfinal= getARGV("-tfinal",0.01);
---
> real Tfinal= getARGV("-tfinal",0.1);
23a24
> Vh ut=u;
30,32c31,42
< solve Pb(u,v) = int2d(Th) ((1/dt)*u*v) 
< + int2d(Th) ( [dx(u),dy(u)]'* [dx(v),dy(v)] - c*u*v) 
< - int2d(Th) ((1/dt)*up*v) - int2d(Th)(f*v) + on(1,2,3,4, u=g);
---
> real cold=c;
> 
> // solve Pb(u,v) = int2d(Th) ((1/dt)*u*v) 
> // + int2d(Th) ( [dx(u),dy(u)]'* [dx(v),dy(v)] - c*u*v) 
> // - int2d(Th) ((1/dt)*up*v) - int2d(Th)(f*v) + on(1,2,3,4, u=g);
> solve Pb(u,v) = 	 
> int2d(Th) (dt*( grad(u)'* grad(v) )) 
> + int2d(Th) (dt*( c*u*v)) 
> -int2d(Th)( ut*v) // TODO check sign lol 
> +int2d(Th)( u*v)
> + on(1,2,3,4, u=g);
> c = sqrt(int2d(Th)(u*u));
38a49
> ut[]=u[];

```

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 4, 2023, 12:23pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/37 "2023-05-04T12:23:24Z")

</div>

First of all, I would like to thank you for the help you have given me. It encourages me to learn more about FreeFem++.  
.  
I am not yet used to using the new version of Freefem  
so I noticed that you have integrated a lot of things in your code.  
My first concern is to know the usefulness of  
include “macro\_ddm.idp”  
and why define the parameters dt and Tfinal in this form.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 4, 2023, 12:32pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/38 "2023-05-04T12:32:22Z")

</div>

How does the output compare to your expectations?

---

<div class="post-metadata">

**Author:** ![aweni](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@aweni](https://community.freefem.org/u/aweni)\
**Post date:** [May 4, 2023, 12:35pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/39 "2023-05-04T12:35:22Z")

</div>

I would also like to know how to generate a file containing E(t) and t in order to plot E(t) with the new Freefem ++4.12

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [May 4, 2023, 1:00pm UTC](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394/40 "2023-05-04T13:00:47Z")

</div>

Do the results look right?

[Previous page](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394.md?page=1)

[Next page](https://community.freefem.org/t/mesh-of-a-right-paving-stone-and-a-cube/2394.md?page=3)