Second order dxx, solver crash

General Discussion
1 
          Hi everybody,

I have a difficulty with the second order  dxx().

Take a simple model, u’’’(x)=1 <=> int2d u’’(x)v(x)=int2d 6v(x). The solution is polynomial of degree 3, like x^3 + constant

mesh Th=square(100,100);
fespace Vh(Th, P2);
Vh u, v;

solve equa(u,v,solver=UMFPACK) =
int2d(Th) ( -dxx(u)dx(v) ) - int2d(Th) ( 6.0v ) + on( 4, u= 1.0 )+ on( 2, u= 2.0 ) ;

But it doesn’'t work. The example bilapMorley.edp for the bilaplacien work perfectly. I modify it to have a 1D problem but the solver crash.

Is it a problem of formulation or the fespace is not correct ?

            Sincerely
2

First, You use 2D finite element for 1d problem,
secondly, you use P2 finite element to compute third derivative
third: you vationnal form is crazy if v = 1 then you have 0 = 6 => system not invertible .

A code in 1d (wrong ) the formulation is wrong pb of Boundary condition
you have 3 constants unknown and you set only 2 …

load "msh3" 
meshL Th=segment(100);
fespace Vh(Th, P2); 
Vh u, v ,xx=x, xxx=x*x*x;

solve equa(u,v,solver=UMFPACK) =
int1d(Th) ( -dxx(u)*dx(v) + 1e-5*dx(u)*dx(v)) - int1d(Th) ( 6*(v) ) + on(1,2,u=xxx);

plot([xx[],u[]],[xx[],xxx[]],wait=1);
3 
  Hi,

I have some setbacks. Your code does’nt work on my machine, is it abnormal or voluntary ?
The error code is :

Build Nodes/DF on mesh : n.v. 101, n. elmt. 100, n b. elmt. 2
nb of Nodes 201 nb of DoF 201 DFon=1100
– FESpace: Nb of Nodes 201 Nb of DoF 201
current line = 11
Assertion fail : (0)
line :245, in file …/femlib/P012_3dCurve.cpp
Assertion fail : (0)
line :245, in file …/femlib/P012_3dCurve.cpp
err code 6 , mpirank 0

I have ever tried to reinstall 2 times from the binaries, version 4.6 and to compile from the source , version 4.10.
Before that, i used freefem 4.4 which don’t use 1 D element, it’s why my code use 2 D element.

I’m much more used to finite difference method , sorry for my dumb question. I don’t understand why P2 finite element is bad in this case ? Because I need dxx so second order ?

3 constants is needed, now i see that, x³+1 is not the only solution without enough constraints.

About the third remark, did you want to say that we have a conflict between int1d(Th)(6*v) and on(2,u=xxx ) (xxx=0) if v=1 at x=0 ?

But applying the integration by parts :

y’’ =6 => y’‘*v =6 *v <=> int2d -y’‘*v’ + int1d y’‘v = int2d 6v <=> int2d y’*v’’ - int1d y’*v’ + int1d y’'v = int2d 6v
I suppose that a factor is missing ? I tried several patchs but nothing work.

      Thank for your reply
4

You have to use the last FreeFem version 4.11 no version 4.6 (Thu Apr 2 14:11:06 2020 +0200)

5

Now it work with the version 4.11 from binary packages, (the last test was under 4.10 ).

I tried to add : +int0d(Th,1) (dxx(u)*v) for dxx(u)=0 at x=0,
or +int0d(Th,2) (dxx(u)v) - int0d(Th,2) ( 6v ) for dxx(u)=6 at x=1
but nothing work.