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 <[email protected]> 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 ([email protected]) 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*
# -*- coding: utf-8 -*-
###########################################################################
################ Yassine CutFEM ##########################
###########################################################################
import getfem as gf
import numpy as np
import matplotlib.pyplot as plt
## Parameters
NX = 8 # Mesh parameter.
Nb_NX = 4
Tab_NX = []
Tab_E_L2 = []
Tab_E_H1 = []
for i in range(Nb_NX):
NX=NX*2
print("The value of NX is")
print(NX)
Id_Omega_h = 1
Id_Gamma_h = 2
Id_Edges = 3
# Parameter du cercle domaine de calcul
x0 = 0.
y0 = 0.
r = 0.34
Coef_gamma = 5.0*NX
# Create a simple cartesian mesh
m = gf.Mesh('triangles grid', np.arange(-.5,.5+1./NX,1./NX),
np.arange(-.5,.5+1./NX,1./NX))
# Create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
mfu = gf.MeshFem(m, 1)
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()
mls.cut_mesh().export_to_pos('MLS.pos')
# assign the P1 fem to all convexes of the both MeshFem
mfu.set_fem(gf.Fem('FEM_PK(2,1)'))
# Integration method used
mim = gf.MeshIm('levelset',mls,'inside', gf.Integ('IM_TRIANGLE(6)'))
mim_bound = gf.MeshIm('levelset',mls,'boundary', gf.Integ('IM_TRIANGLE(6)'))
mim.set_integ(6)
#mfu.dof_from_im(mim)
# Reperage des domaines
Tab_CNV_Ind_Omega = mim.convex_index()
m.set_region(Id_Omega_h, m.faces_from_cvid(Tab_CNV_Ind_Omega))
#print("Omega_h")
#print(m.region(Id_Omega_h))
# Set the region Gamma_h
m.set_region(Id_Gamma_h, m.outer_faces(Tab_CNV_Ind_Omega))
#print("Gamma_h")
#print(m.region(Id_Gamma_h))
# Set the list of internal edges for the Ghost penalty
Tab_band_Gamma = mim_bound.convex_index()
Tab_all_ed_band = m.faces_from_cvid(Tab_band_Gamma)
m.set_region(Id_Edges,Tab_all_ed_band)
m.set_region(6,m.outer_faces(Tab_CNV_Ind_Omega))
m.region_subtract(Id_Edges, 6)
#print("the liste of inter edges in the band")
#print(m.region(Id_Edges))
# Model
md = gf.Model('real')
# Main unknown
md.add_fem_variable('u', mfu)
# Interpolate the exact solution (Assuming mfu is a Lagrange fem)
Uexpr = ('(np.sqrt(np.square(x-%f)+np.square(y-%f))**3-%f**3)/9.0'%(x0,y0,r));
Uex = mfu.eval(Uexpr,globals(),locals())
# Laplacian term on u
md.add_Laplacian_brick(mim, 'u')
# Volumic source term
Fexpr = ('-np.sqrt(np.square(x-%f)+np.square(y-%f))'%(x0,y0));
F = mfu.eval(Fexpr,globals(),locals())
md.add_initialized_fem_data('f', mfu, F)
md.add_source_term_brick(mim, 'u', 'f')
# Add the first Boundary condition on Gamma_h
Var_Normal_Grad_u = "Grad_u.Normal"
Test_of_u = "Test_u"
md.add_linear_term(mim_bound, "-({S})*({K})".format(S=Var_Normal_Grad_u, K=Test_of_u), -1)
#Add the second Boudary term on Gamma
var_u = "u"
Normal_test_of_u = "Grad_Test_u.Normal"
md.add_linear_term(mim_bound, "-({L})*({G})".format(L=var_u, G=Normal_test_of_u), -1)
# Add the 3th condition with penalization
#g = mfu.eval('0')
#md.add_initialized_fem_data('g', mfu, g)
md.add_initialized_data('g', [0.0])
md.add_Dirichlet_condition_with_penalization(mim_bound, 'u', Coef_gamma, -1, 'g', mfu)
# Add the 4th Boundary condition Ghost Penalty
md.add_initialized_data('Sigma_h', [0.1/NX])
jump = "((Grad_u-Interpolate(Grad_u,neighbour_elt)).Normal)"
test_jump = "((Grad_Test_u-Interpolate(Grad_Test_u,neighbour_elt)).Normal)"
md.add_linear_term(mim_bound, "Sigma_h*({M}).({N})".format(M=jump, N=test_jump), Id_Edges)
#md.add_Dirichlet_condition_with_penalization(mim_bound,
#md.add_linear_term(mim_bound, 'Sigma_h*(Xfem_plus(Grad_u)-Xfem_minus(Grad_u)).(Xfem_plus(Grad_Test_u)-Xfem_minus(Grad_Test_u))', Id_Edges)
# Assembly of the linear system and solve.
md.solve()
# Main unknown
Uap = md.variable('u')
# Slice pour la restriction de la solution et le calcul de l'erreur
sl = gf.Slice(('comp',('ball', +1, [x0, y0], r)), m, 5)
# Restriction of the solution to the physical domain
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)
# The approximate solution on the fictititious domain
mfu.export_to_pos('Solution.pos', Uap,'Approximate solution')
mfu.export_to_pos('ExSolution.pos', Uex,'Exact solution')
mfu.export_to_pos('Poisson-Exact-App.pos', Uex,'Exact solution',
Uap,'Computed solution')
L2error = gf.compute(mfu, Uap-Uex, 'L2 norm', mim)
H1error = gf.compute(mfu, Uap-Uex, 'H1 norm', mim)
print('Error in L2 norm : ', L2error)
print('Error in H1 norm : ', H1error)
# Sauvegarde des résultats pour le calcul de l'ordre de convergence
Tab_NX.append(np.log(1./NX))
Tab_E_L2.append(np.log(L2error))
Tab_E_H1.append(np.log(H1error))
# plot the convergence rate
plt.plot(Tab_NX, Tab_E_H1, 'v-', Tab_NX, Tab_E_L2, 'o-')
plt.legend(["H1 error", "L2 error"])
plt.axes()
plt.xlabel("Log of the meshsize")
plt.ylabel('Log of the error')
plt.show()
# Calcul d'ordre de convergence pour les deux normes
Order_L2 =[]
Order_H1 =[]
for i in range(1,Nb_NX):
OL2 = (Tab_E_L2[i]-Tab_E_L2[i-1])/(Tab_NX[i]-Tab_NX[i-1])
OH1 = (Tab_E_H1[i]- Tab_E_H1[i-1])/(Tab_NX[i]-Tab_NX[i-1])
Order_L2.append(OL2)
Order_H1.append(OH1)
print("L'ordre de convergence selon la norme L2")
print(Order_L2)
print("L'ordre de convergence selon la norme H1")
print(Order_H1)
