#===============================================================================================================================
# Name : Résolution EF de l'équation de Laplace-Beltrami -\triangle u = f sur la frontière d'un cube
# Author : Michael Ndjinga
-# Copyright : CEA Saclay 2020
+# Copyright : CEA Saclay 2021
# Description : Utilisation de la méthode des éléménts finis P1 avec champs u et f discrétisés aux noeuds d'un maillage triangulaire
# Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
# Résolution d'un système linéaire à matrice singulière : les vecteurs constants sont dans le noyau
-# Comparaison de la solution numérique avec la solution exacte définie face par face :
-# u= sin(2*pi*x)*sin(2*pi*y) (haut et bas)
-# u=-sin(2*pi*x)*sin(2*pi*z) (gauche et droite)
-# u= sin(2*pi*y)*sin(2*pi*z) (avant et arrière)
+# Comparaison de la solution numérique avec la solution exacte définie face par face : u(x,y,z)= cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
#================================================================================================================================
import cdmath
-from math import sin, pi
+from math import cos, pi
import numpy as np
import PV_routines
import VTK_routines
y = Ni.y()
z = Ni.z()
- if abs(x)<eps or abs(x-1)<eps: #AVANT et ARRIERE
- my_RHSfield[i]= 8*pi*pi*sin(2*pi*y)*sin(2*pi*z)
- elif abs(y)<eps or abs(y-1)<eps: #GAUCHE et DROITE
- my_RHSfield[i]=-8*pi*pi*sin(2*pi*x)*sin(2*pi*z) #Minus sign to have a smooth function of the curvilinear variable
- elif abs(z)<eps or abs(z-1)<eps: #HAUT et BAS
- my_RHSfield[i]= 8*pi*pi*sin(2*pi*x)*sin(2*pi*y)
- else:
- raise ValueError("Domain should be the unit cube skin with 6 faces")
+ my_RHSfield[i]= 8*pi*pi*cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
if my_mesh.isBorderNode(i): # Détection des noeuds frontière
raise ValueError("Mesh should not contain borders")
# Storing of numerical errors and mesh sizes
for filename in meshList:
error_tab[i], mesh_size_tab[i], min_sol_num, max_sol_num, time_tab[i] =FiniteElements3DPoissonCubeSkin.solve(mesh_path+filename, resolution,meshType,testColor)
- assert min_sol_num>-1.1
- assert max_sol_num<1.1
+ assert min_sol_num>-5.
+ assert max_sol_num<6.
error_tab[i]=log10(error_tab[i])
time_tab[i]=log10(time_tab[i])
- with open('./FiniteElementsOnCubeSkinPoisson_PlotOnSortedLines'+meshType+str(mesh_size_tab[i])+'.csv') as f:
- lines = f.readlines()
- y = [float(line.split(",")[0]) for line in lines[1:]]
- x = [float(line.split(",")[1]) for line in lines[1:]]
+ # with open('./FiniteElementsOnCubeSkinPoisson_PlotOnSortedLines'+meshType+str(mesh_size_tab[i])+'.csv') as f:
+ # lines = f.readlines()
+ # y = [float(line.split(",")[0]) for line in lines[1:]]
+ # x = [float(line.split(",")[1]) for line in lines[1:]]
- plt.plot(x, y, label= str(mesh_size_tab[i]) + ' nodes')
+ # plt.plot(x, y, label= str(mesh_size_tab[i]) + ' nodes')
mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
i=i+1
b=(-a2*b1+a1*b2)/det
print( "FE on 3D cube skin triangle mesh : scheme order is ", -a)
- assert abs(a+0.816)<0.1
+ assert abs(a+1.915)<0.001
# Plot of convergence curves
plt.close()
#===============================================================================================================================
# Name : Résolution EF de l'équation de Laplace-Beltrami -\triangle u = f sur la frontière d'un cube
# Author : Michael Ndjinga
-# Copyright : CEA Saclay 2020
+# Copyright : CEA Saclay 2021
# Description : Utilisation de la méthode des éléménts finis P1 avec champs u et f discrétisés aux noeuds d'un maillage triangulaire
# Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
# Résolution d'un système linéaire à matrice singulière : les vecteurs constants sont dans le noyau
-# Comparaison de la solution numérique avec la solution exacte définie face par face :
-# u=sin(2*pi*x)*sin(2*pi*y) (haut et bas)
-# u=sin(2*pi*x)*sin(2*pi*z) (gauche et droite)
-# u=sin(2*pi*y)*sin(2*pi*z) (avant et arrière)
+# Comparaison de la solution numérique avec la solution exacte définie face par face : u(x,y,z)= cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
#================================================================================================================================
import cdmath
import time, json
-from math import sin, pi
+from math import cos, pi
import PV_routines
import VTK_routines
import paraview.simple as pvs
y = Ni.y()
z = Ni.z()
- if abs(x)<eps or abs(x-1)<eps: #AVANT et ARRIERE
- my_RHSfield[i]= 8*pi*pi*sin(2*pi*y)*sin(2*pi*z)
- elif abs(y)<eps or abs(y-1)<eps: #GAUCHE et DROITE
- my_RHSfield[i]=-8*pi*pi*sin(2*pi*x)*sin(2*pi*z) #Minus sign to have a smooth function of the curvilinear variable
- elif abs(z)<eps or abs(z-1)<eps: #HAUT et BAS
- my_RHSfield[i]= 8*pi*pi*sin(2*pi*x)*sin(2*pi*y)
- else:
- raise ValueError("Domain should be the unit cube skin with 6 faces")
+ my_RHSfield[i]= 8*pi*pi*cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
if my_mesh.isBorderNode(i): # Détection des noeuds frontière
raise ValueError("Mesh should not contain borders")
finiteElementsOnCubeSkin_0vtuDisplay = pvs.Show(finiteElementsOnCubeSkin_0vtu, renderView1)
pvs.ColorBy(finiteElementsOnCubeSkin_0vtuDisplay, ('POINTS', 'ResultField'))
slice1Display = pvs.Show(slice1, renderView1)
- pvs.SaveScreenshot("./FiniteElementsOnCubeSkinPoisson"+"_Slice_"+meshType+str(nbNodes)+'.png', magnification=1, quality=100, view=renderView1)
+ pvs.SaveScreenshot("./FiniteElementsOnCubeSkinPoisson"+"_Slice_"+meshType+str(nbNodes)+".png", magnification=1, quality=100, view=renderView1)
plotOnSortedLines1 = pvs.PlotOnSortedLines(Input=slice1)
pvs.SaveData('./FiniteElementsOnCubeSkinPoisson_PlotOnSortedLines'+meshType+str(nbNodes)+'.csv', proxy=plotOnSortedLines1)
