Re: [Getfem-users] Calculate the error in the physical domain

[Yassine ZAIM]

Dear Yves Renard,

Thank you so much for your several appreciated helps. 

I will return to your previous remark about the Ghost penalty, you are right, I confirm that the normal derivative is sufficient for the P2 and P3 methods in addition to the alternative method proposed in your paper.

I tried different values of the parameters (penalty coefficients) bigger and smaller. But unfortunately, I don't get the optimal rate of convergence. The error is approximately the same for the different values of the parameters, it takes the minimum values near the boundary and the maximum values in the center of the circle. 

I did also two other tests with big numbers of subdivisions NX=600 and NX=1200 and I get in this case the following rate of convergence :

Order with respect to H1 norm = 1.4967550587246639

Order with respect to L2 norm = 1.997064049265697

But it is not the case for the small values of NX. 

Best regards. 

Le mer. 2 janv. 2019 à 13:41, Yves Renard <address@hidden> a écrit :

Dear Yassine,

Ok, it seems to me that your error computation is correct and should give you the correct convergence rate.
Did you try to adapt your penalty coefficient (it can lock if it is too high) and try to see where the error is concentrated ?


Best regards,

Yves.



----- Original Message -----
From: "Yassine ZAIM" <address@hidden>
To: "yves renard" <address@hidden>
Cc: "getfem-users" <address@hidden>
Sent: Monday, December 31, 2018 10:11:59 AM
Subject: Re: [Getfem-users] Calculate the error in the physical domain

Dear Yves Renard,
Thank you so much for your answers. The alternative to the Ghost penalty is
a good idea to solve with P_k FEMs. For the remark about the slice, I don't
use it to calculate the error, but just to export the solution in the
physical domain. To calculate the error I use exactly the integration
method with  MeshLevelSet as follow:
ls = gf.LevelSet(m,2,'(x-{X0})*(x-{X0}) + (y-{Y0})*(y-{Y0}) -
{R}*{R}'.format(X0=x0, Y0=y0, R=r))
mls = gf.MeshLevelSet(m)
mls.add(ls)
mls.adapt()
mim = gf.MeshIm('levelset',mls,'INSIDE', gf.Integ('IM_TRIANGLE(6)'))
mim.set_integ(6)
#After creating the model and solve it, I calculate the error as follows:#
    L2error = gf.compute(mfu, Uap-Uex, 'L2 norm', mim)
    H1error = gf.compute(mfu, Uap-Uex, 'H1 norm', mim)
where mfu is given by:
mfu   = gf.MeshFem(m, 1)
mfu.set_fem(gf.Fem('FEM_PK(2,1)'))
But I still not able to get the optimal rate of convergence (I get
approximately the same error for NX=16, 32, 64,etc. But if I take NX=400
for example, I get an approximate solution very similar to the exact
solution with a very small error where the maximum value=2.1e-06).
Thank you for the help.

Le jeu. 27 déc. 2018 à 09:54, Yves Renard <address@hidden> a
écrit :

>
> Dear Yassine,
>
> Concerning the Ghost penalty, I am not sure that penalization of normal
> derivative is not sufficient even for P2 or P3 method. I think it should be
> sufficient.
>
>
> Note also that an alternative to Ghost penalty is to directly interpolate
> the normal derivative in a neighbour element for elements with a small
> intersection with the domain as we proposed in the paper
>
> J. Haslinger, Y. Renard. A new fictitious domain approach inspired by the
> extended finite element method. Siam J. on Numer. Anal., 47(2):1474--1499,
> 2009.
> (http://math.univ-lyon1.fr/~renard/papers/2009_fictif.pdf)
>
> Best regards,
>
> Yves
>
>
> ----- Original Message -----
> From: "yves renard" <address@hidden>
> To: "Yassine ZAIM" <address@hidden>
> Cc: "getfem-users" <address@hidden>
> Sent: Wednesday, December 26, 2018 10:02:11 PM
> Subject: Re: [Getfem-users] Calculate the error in the physical domain
>
> Dear Yassine,
>
> Once again, the slice operations are only for graphical post-processing.
> You cannot compute an error with them. You have to use cut integration
> methods instead.
>
> Best regards,
>
> Yves
>
> ----- Original Message -----
> From: "Yassine ZAIM" <address@hidden>
> To: "yves renard" <address@hidden>
> Cc: "getfem-users" <address@hidden>
> Sent: Wednesday, December 19, 2018 5:00:01 PM
> Subject: Re: [Getfem-users] Calculate the error in the physical domain
>
> Dear Yves Renard and getfem++ users,
> Thank you so much for your answer. I am inspired by the
> "demo_fictitious_domains.py"
> example to achieve my test and the proposed method is what I did at the
> beginning. But I get the same value of the error for different values of NX
> and I can't get the optimal order. For this, I am saying maybe the outside
> domain of Omega which creates the problem, for this I asked how to restrict
> the computation of error to Omega. I enclosed my code if you can see it
> quickly, maybe you will find out the mistake.
> I have another question about the Ghost penalty. For the P_2 element, I
> guess that I have to work with the hessian "hess_u" and "hess_Test_u" But
> how we can implement it for P_3, P_4,... spaces.
> Thank you so much for your help.
>
> Le mer. 19 déc. 2018 à 14:39, Yves Renard <address@hidden> a
> écrit :
>
> >
> > Dear Yassine,
> >
> > The slice operation is only for graphical post-processing. You cannot use
> > a sliced solution into computations because the obtained vector represent
> > the value of the solution on a cut mesh representing the computed
> > intersection. The obtained vector is no longer a dof vector on the
> meshfem
> > mfu. The standard way to compute this error is to produce a cut
> integration
> > method, and just compute the L2 or H1 norms in a standard way with it.
> You
> > can produce cut integration method with the MeshLevelSet and
> > MeshIM('levelset') objects. An example of use in the python test program
> > "demo_fictitious_domains.py"
> >
> > Best regards,
> >
> > Yves
> >
> >
> > Le 18/12/2018 à 18:54, Yassine ZAIM a écrit :
> >
> > Dear getfem++ users,
> > I am trying to learn how to implement the fictitious domain method for
> the
> > simple problem of Poisson. My domain of interest Omega is a circle of
> > center (*x0=0,y0=0*) and radius *r=0.35*. I made the resolution in the
> > fictitious domain (the square [-0.5, 0.5]^2) and after that, I used the
> > slice to restrict my solution to the physical domain.
> >     sl = gf.Slice(('comp',('ball', +1, [x0, y0], r)), m, 5)
> >     sl.export_to_vtk('App_Solution.vtk', 'ascii', mfu, Uap)
> >     sl.export_to_vtk('Exc_Solution.vtk', 'ascii', mfu, Uex)
> >     sl.export_to_vtk('Error-Ex-App.vtk', 'ascii', mfu, Uap-Uex)
> > When I see the approximate and the exact solution in addition to the
> error
> > I can say that I get a good approximate solution. But I am trying to get
> > the optimal order of convergence (h and h^2) like in the papers "J.
> > Haslinger and Y. Renard" or "E. Burman and P. Hansbo". So I interpolate
> the
> > exact and approximate solution in the slice (domain of interest).
> >     Ue = gf.compute(mfu, Uex, 'interpolate on', sl)
> >     U = gf.compute(mfu, Uap, 'interpolate on', sl)
> > And after that I tried to calculate the error for different values of
> > NX=[16,32,..]
> >     L2error = gf.compute(mfu, U-Ue, 'L2 norm', mim)
> >     H1error = gf.compute(mfu, U-Ue, 'H1 norm', mim)
> > such that :
> > mim = gf.MeshIm('levelset',mls,'inside', gf.Integ('IM_TRIANGLE(6)'))
> > and  mfu.set_fem(gf.Fem('FEM_PK(2,1)'))
> > But I get the following error in the "compute" function (of the errors) :
> > return getfem('compute', mf, U, what, *args)
> > RuntimeError: (Getfem::InterfaceError) -- The trailing dimension of
> > argument 2 (an array of size 4670) has 4670 elements, 289 were expected.
> > My question is how to calculate the error just in the interesting domain
> > (physical domain, slice sl in my case).
> > Thank you in advance.
> >
> >
> > --
> > *ZAIM Yassine *
> > *PhD in Applied Mathematics*
> >
> >
> > --
> >
> >   Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
> >   INSA-Lyon
> >   20, rue Albert Einstein
> >   69621 Villeurbanne Cedex, FRANCE
> >   http://math.univ-lyon1.fr/~renard
> >
> > ---------
> >
> >
>
> --
> *ZAIM Yassine *
> *PhD in Applied Mathematics*
>


--
*ZAIM Yassine *
*PhD in Applied Mathematics*

--

ZAIM Yassine 

PhD in Applied Mathematics